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ABSTRACT 

This  investigation  develops  an  axisymmetric  heat  transfer- 
combustion  model  of  a  porous  medium  within  a  circular  cylin- 
der.  System  flow  is  governed  by  Darcy's  law.   Carbon  and 
air  properties  are  treated  as  variables  of  temperature.   A 
combined  continuity-Darcy  equation,  an  oxygen  mass  balance 
equation,  and  energy  balance  equations  (one  each)  for  air 
and  carbon,  describe  the  conservation  laws  of  the  system. 
Transport  mechanisms  for  oxygen  mass  transfer  are  molecular 
diffusion  and  convective  transport,  and  an  oxygen  consumption 
term  to  account  for  combustion  is  included.   Heat  transfer 
mechanisms  included  in  the  model  are  conduction  and  convec- 
tion.  Radiation  is  accounted  for  at  applicable  boundaries 
only.   Nonvolatile  combustion  is  accounted  for  in  the  carbon 
energy  and  oxygen  mass  balance  equations  as  a  heat  generation 
term  of  Arrhenius  type.   The  numerical  solution  of  four 
coupled,  nonlinear,  transient  partial  differential  field 
equations  is  accomplished  using  the  Galerkin  formulation  of 
the  Finite  Element  Method.   The  effect  of  porosity  on  system 
behavior  is  examined. 
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I.   INTRODUCTION 

A.   PRIOR  INVESTIGATIONS 

The  problem  of  characterizing  physical  systems  involving 
combustion  and  quantifying  the  accompanying  thermal  response 
has  received  attention  in  recent  years.   The  level  of  diffi- 
culty encountered  in  a  problem  of  this  type  precluded 
analytical  as  well  as  numerical  solutions.   The  difficulty 
of  the  problem  lies  in  the  number  of  disciplines  that  encom- 
pass it.   The  problem  involves  the  kinetics  of  combustion, 
heat  and  mass  transfer  mechanisms  and  fluid  flow.   Numerical 
solutions  with  increased  efficiency  in  computation  are  now 
possible  with  high  speed  com.puters. 

The  combustion  problem  has  applications  in  the  areas  of 
forest  fire  control,  energy  conservation,  underground 
(nuclear)  fuel  storage  and  waste  disposal,  and  natural  gas 
fire  control.   Major  contributions  may  be  made  in  the  area 
of  energy  production  and  conservation  by  the  analysis  of 
heat  generation  and  ignition  characteristics  of  combustible 
materials . 

The  characterization  of  combustion  and  heat  transfer  in 
porous  media  has  been  of  particular  interest.   A  brief  survey 
of  some  of  the  investigations  in  the  area  of  combustion  and 
heat  transfer  in  porous  media  follows.   The  intent  is  to 
acquaint  the  reader  with  work  in  the  field  that  offered 
insight  into  the  formulation  of  this  investigation. 
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Investigations  by  Kordylewski  on  the  influence  of  aero- 
dynamics on  the  critical  parameters  of  thermal  ignition 
[Ref.  1],  examined  homogeneous  combustion  in  a  two-dimensional 
cylindrical  reactor.   The  transient  analysis  involved 
Arrhenius  combustion  of  a  solid  fuel.   Heat  transfer  mechanisms 
considered  were  convection  and  conduction.   Because  the  inves- 
tigation dealt  with  thermal  ignition  theory,  reactant  con- 
sumption was  omitted.   The  flow  field  was  assumed  to  be 
steady  prior  to  ignition  and  constant  fluid  properties  were 
assumed. 

Safety  dictates  the  assessment  of  the  structural  integrity 
of  a  building  after  a  severe  fire.   In  order  to  predict 
stresses  due  to  fire  in  buildings,  the  thermal  response  must 
be  known.   Sahota  and  Pagni  [Ref.  2],  formulated  a  transient 
solution  for  two-phase,  two-component  flow  in  one-dimensional 
porous  concrete  structures.   The  mechanisms  considered  in- 
cluded:  heat  conduction,  molecular  diffusion  of  gaseous 
components,  and  pressure-driven  convective  flow  subject  to 
Darcy's  law. 

A  sudden  reduction  in  feed  temperature  in  a  packed-bed 
reactor  leads  to  the  transient  temperature  rise  known  as 
"wrong-way  behavior."  This  behavior  was  investigated  by 
Mehta,  Sams,  and  Luss  [Ref.  3].  The  work  identifies  the 
important  rate  processes  and  parameters  which  cause  the 
behavior,  and  generates  an  expression  for  predicting  the 
magnitude  of  the  maximum  transient  temperature. 
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Ignition  parameters  in  porous  solid  fuels  have  been 
analyzed  by  Kim  and  Chung  [Ref.  4].   They  investigated 
three  porous  solid  fuel  geometries  (a  semi-infinite  slab,  an 
infinitely  long  circular  cylinder  and  a  sphere)  with  con- 
stant energy  and  gaseous  oxidant  fluxes  at  the  fuel  surface. 
Laplace  transformation  of  the  nondimensionalized  oxidant  mass 
equation  and  fuel  energy  equation  allowed  for  asymptotic 
solution  of  a  nondimensionalized  transient  temperature 
equation  which  is  valid  in  the  neighborhood  of  the  fuel 
surface.   Observations  included  shorter  ignition  times  for 
spheres  compared  to  slabs.   Times  for  ignition  increased  with 
fuel  size  and  approached  values  of  semi-infinite  slabs 
asymptotically.   Other  effects  of  size  and  geometry  of  porous 
solid  fuels  on  ignition  parameters  are  presented. 

Saatdjian  and  Caltagirone  [Ref.  5]  investigated  a  transient 
two-dimensional  combustion  model  with  natural  convection. 
A  porous  medium  undergoing  exothermic  combustion  was  saturated 
by  a  gas  and  bounded  by  two  impermeable  planes.   Permeability, 
fluid  viscosity,  thermal  conductivity  of  the  porous  medium, 
and  the  heat  of  reaction  were  assumed  to  be  constants. 

Home  and  O'Sullivan  [Ref.  6],  Hickox  [Ref.  7]  and  Chan 
and  Banerjee  [Ref.  8],  investigated  the  natural  convection 
phenomenon  in  porous  media.   Exothermic  reactions  were  not 
addressed  in  these  investigations.   Home  and  O'Sullivan 
considered  the  effects  of  variable  viscosity  on  the  stability 
of  a  porous  layer  in  a  transient  two-dimensional  problem. 
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The  boundary  configuration  for  this  work  was  a  two-dimensional 
region  with  a  line  of  vertical  symmetry,  enclosed  on  both 
sides  by  impermeable,  insulated  sidewalls,  and  heated  along 
half  of  the  base.   The  convective  flow  of  fluid  through 
porous  media  heated  from  below  has  applications  in  the  study 
of  behavior  of  geothermal  systems.   Hickox  studied  convection 
as  an  application  to  sub-seabed  disposal  of  nuclear  waste. 
The  problem  was  a  two-dimensional  transient  analysis  of  free 
convection  produced  by  a  concentrated  heat  source  (implanted 
container  of  waste  material)  in  the  subsurface  sedimentary 
layer  of  a  seabed.   The  porous  medium  was  assumed  to  be  rigid, 
homogeneous  and  isotropic.   Density  changes  of  the  fluid  were 
accounted  for  only  in  the  buoyancy  term  of  Darcy's  law. 
Permeability,  viscosity  and  thermal  conductivity  were  assumed 
to  be  constant.   Chan  and  Banerjee  conducted  a  transient 
three-dimensional  analysis  of  natural  convection  in  porous 
media.   In  addition  to  convective  heat  transfer,  conduction 
was  also  considered  between  solid  spherical  particles  that 
comprise  the  porous  medium.   The  porous  medium  was  considered 
homogeneous  and  isotropic  in  its  physical  properties  including 
permeability  and  thermal  conductivity,  both  of  which  were 
assumed  to  be  constant  with  temperature.   Fluid  density  was 
considered  to  be  constant  except  in  the  buoyancy  term  of 
Darcy's  law. 

In  1982,  Vatikiotis  [Ref.  9]  considered  a  transient  one- 
dimensional  heat  transfer  and  combustion  model  of  a  porous 
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medium.   The  heat  transfer  mechanisms  included  were  conduc- 
tion, convection  and  radiation.   One  problem  considered  was 
a  porous  mat  subject  to  Arrhenius  combustion  in  which  all 
properties  were  temperature  dependent.   The  capillary  serial 
(permeability)  model  introduced  by  Scheidegger  [Ref.  10]  was 
employed. 

The  current  investigation  considers  the  effects  of 
porosity  on  system  behavior  for  a  two-dimensional  (axi- 
symmetric)  model  of  combustion  in  a  porous  medium.   A  great 
part  of  the  following  analysis  may  be  found  in  Vatikiotis 
[Ref.  9] .   The  pertinent  parts  are  repeated  here  for  conven- 
ience to  the  reader.   Deviations  are  cited  and  are  for  the 
most  part  due  to  the  axisymmetric  two-dimensional  aspect  of 
the  present  model  vice  the  previous  one-dimensional  model. 
The  storage  scheme  employed  by  the  numerical  model  warrants 
the  reader's  attention.   The  savings  in  computer  storage 
realized  by  utilizing  the  method  of  Franke  and  Salinas 
[Ref.  29],  is  substantial  and  is  addressed  in  the  Numerical 
Section  of  this  work.   The  nonlinear  combustion  term  is 
treated  as . a  bilinear  spatial  operator  of  carbon  temperature 
and  oxygen  concentration,  in  contrast  to  the  Vatikiotis  treat- 
ment of  the  term  as  an  excitation  vector.   The  idea  behind 
the  present  treatment  was  to  capture  the  effects  of  the  non- 
linear combustion  term  on  both  of  these  dependent  variables. 
It  was  felt  that  this  would  alleviate  some  numerical  diffi- 
culties.  A  numerical  model  is  formulated  and  results  are 
presented. 
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II.   PROBLEM  DESCRIPTION 

The  problem  under  investigation  is  that  of  combustion 
of  a  porous  fuel  (carbon)  imbedded  in  a  cylindrical  container. 
Pressure  gradients  induce  convective  currents  through  the 
medium.   The  air  flow  produces  two  opposing  effects: 
(1)  internal  convective  heat  transfer,  and  (2)  a  supply  of 
oxygen  to  promote  heat  generation  through  combustion. 
Extinguishment  or  sustained  combustion  of  the  porous  medium 
depend  on  the  interaction  of  these  effects.   If  heat  transfer 
dominates,  the  combustion  will  proceed  to  extinguishment, 
otherwise  combustion  will  continue.   A  mathematical  model  was 
developed  to  provide  an  understanding  of  this  interaction 
and  its  effect  on  thermal  behavior. 

The  mathematical  model  is  formulated  as  follows.   Energy 
balances  on  the  carbon  and  convected  air  provide  heat 
transfer  equations  for  each.   The  heat  transfer  mechanisms 
incorporated  in  the  model  are:   (1)  conduction,  (2)  convection 
and  (3)  radiation  (where  applicable,  at  boundaries).   In 
addition,  nonvolatile  combustion  is  included  in  the  carbon 
heat  transfer  equation  as  a  fractional  order  heat  generation 
term  of  Arrhenius  type. 

The  conservation  of  species  law  is  applied  to  the  oxygen 
molecule  concentration  to  yield  a  third  equation.   The 
resulting  oxygen  mass  transfer  equation  includes  the 
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transfer  mechanism  of:   (1)  molecular  diffusion,  and 

(2)  convective  transport  of  species.   An  oxygen  consumption 

term  due  to  combustion  is  included. 

The  fourth  field  equation  involves  the  system  pressure 
gradient.   The  equation  is  a  combined  Darcy's  law  and  air 
mass  continuity  equation  for  the  system. 

The  conservation  equations  describing  the  system  field 
are  four  transient,  coupled,  nonlinear  partial  differential 
equations.   The  four  equations  are  solved  by  a  two-dimensional 
Galerkin  formulation  of  the  Finite  Element  Method.   The  inte- 
gration scheme  used  for  this  highly  stiff  system  is  one 
presented  by  Gear  and  modified  by  Franke  [Ref.  30].   The 
scheme  is  especially  suited  to  systems  of  implicit  and  stiff 
differential  equations.   The  integration  scheme  is  used  in 
conjunction  with  an  optimal  compact  storage  scheme  proposed 
by  Franke  and  Salinas  [Ref.  29j . 
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III.   THEORY  AND  BACKGROUND 

A.   DESCRIPTION  OF  THE  POROUS  MEDIUM 

In  this  work,  a  porous  medium  is  considered  to  be  a 
solid  containing  interconnecting  pores  that  allow  fluid  to 
permeate  and  flow  through  the  solid.   Either  of  two  classes 
of  porous  media,  consolidated  (solid  and  rigid)  and  uncon- 
solidated (comprised  of  discrete  particles  as  found  in 
granular  beds)  may  be  considered.   Each  class  of  porous 
media  may  have  isotropic  nonhomogeneous  properties. 

The  common  characteristics  of  all  porous  media  are: 
(1)  porosity,  (2)  specific  internal  area,  (3)  pore  diameter, 
and  (4)  tortuosity.   Porosity,  p,  is  defined  as  the  ratio  of 
void  volume  to  total  volume.   The  specific  internal  area,  z, 
is  the  ratio  of  internal  surface  area  to  bulk  volume.   In 
general,  the  distribution  of  pore  size  in  a  porous  medium  is 
random  (nonhomogeneous)  and  dynamic  (subject  to  small 
strains  induced  by  the  transient  pressure  field) .   The  situa- 
tion motivates  the  investigator  to  treat  the  porous  medium 
as  a  continuum  possessing  idealistic  geometrical  properties 
of  porosity  and  pore  diameter.   The  specific  internal  area 
is  obtained  from  the  porosity  model.   Though  many  conventions 
exist  for  the  definition  of  a  conceptual  pore  diameter,  in 
this  study  a  hydraulic  diameter  analogy  is  employed.   The 
tortuosity,  x,  of  a  porous  medium  is  the  ratio  of  the  flow 
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path  to  the  straight  (line)  path  displacement  of  the 
particle.   Permeability,  a  measure  of  hydraulic  conductivity, 
is  a  property  of  the  porous  medium  that  depends  upon  the 
four  characteristics  mentioned  above.   Scheidegger  [Ref.  10] 
presents  methods  used  to  measure  the  properties  of  porous 
media.   Methods  discussed  are  essentially  experimental  in 
nature. 

The  porous  medium  was  modelled  as  shown  in  Figure  3.1. 
In  Figure  3.1,  D  is  the  particle  center-to-center  distance, 
and  d  is  the  particle  diameter.   From  the  idealized  geometry, 
the  porosity  for  spherical  particles  is. 


P   =   1  -  f(§)'  (3.1 


The  pore  diameter  is  obtained  from  an  expression  proposed 
by  Carman  [Ref.  11], 


6   =   ^  (3.2) 


Equation  3.2  is  analogous  to  a  more  familiar  form  of  mean 
hydraulic  diameter,  4v/A,  where  v  is  the  void  volume  and  A 
is  the  wetted  surface  area.   The  specific  internal  area,  Z, 
based  on  the  idealized  geometry  of  Figure  3.1  may  be  expressed 
as , 


1   h2 
Z   =   yTT^  (3.3) 

^   D 
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for  spherical  particles.   Equation  3.3  is  sometimes  used 
and  assumes  one-half  of  the  total  internal  surface  area 
is  effective  for  convective  heat  transfer.   The  fractional 
amount  of  total  area  is  an  estimate  based  on  Fontenot's  [Ref.  37' 
experimental  results  and  does  not  generally  apply  to  porous 
media  [Ref.  9].   The  Kozeny  relations,  alternate  expressions 
of  the  specific  internal  area,  are  discussed  by  Scheidegger 
[Ref.  10].   Advantages  of  the  Kozeny  relations  are  fair 
agreement  with  experimental  values  and  calculations  that  are 
independent  of  particle  shape.   A  disadvantage  is  the 
failure  of  the  relations  to  predict  accurate  values  of  z  at 
high  values  of  porosity  [Ref.  9].   For  the  geometric 
configuration  of  Figure  3.1,  the  tortuosity  depends  on  the 
ratio  d/D.   Carman  [Ref.  11]  presents  a  table  of  measured 
tortuosity  factors  of  various  materials  and  geometries,  and 
points  out  differences  between  analytical  determinations  of 
tortuosity.   In  this  study,  his  recommended  value  of  1.4  is 
used.   Particle  size  decreases  with  consumption.   As  a 
result,  thermophysical  properties  which  depend  on  particle 
diameter,  as  well  as  temperature,  are  functions  of  time- and 
space.   The  numerical  model  presented  assumes  matrix 
rigidity  as  particle  diameter  decreases.   This  assumption 
gives  rise  to  an  increase  in  porosity  during  combustion. 
Although  Scheidegger  [Ref.  10]  in  his  discussion  of  the 
packing  theory  of  spheres,  reports  .875  porosity  as  the 
threshold  for  stability  in  a  porous  matrix.  Carman  [Ref.  11] 
reports  on  investigations  performed  on  porosities  as  high 
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as  0.99.   Changes  in  particle  diameter  are  incorporated  in 
the  model  and  are  discussed  in  Section  III.D  with  carbon 
combustion. 

B.   DARCY'S  LAW  AND  PORE  VELOCITY 

The  Reynolds  number  for  porous  media  is  defined  by 


p  s  d 
Re   =  — (3.4) 


where  s  is  the  local  pore  velocity,  p  is  the  mass  density  of 
air,  and  y  is  the  dynamic  viscosity.   The  magnitude  of  the 
Reynolds  number  indicates  whether  fluid  motion  is  dominated 
by  molecular,  viscous,  or  inertial  effects.   Most  investi- 
gations of  flow  through  porous  media  indicate  flow  regimes 
of  viscous  and  inertially  dominated  flows.   The  Navier- 
Stokes  equations  apply  to  fluid  motion  possessing  such 
Reynolds  numbers.   Because  of  the  geometry  involved  in  a 
consolidated  (rigid)  porous  medium  and  the  no-slip  boundary 
condition  (i.e.,  s  =  0  at  a  solid-fluid  interface),  solution 
of  the  Navier-Stokes  equations  is  difficult  for  a  porous 
medium.   Scheidegger  [Ref.  10]  points  out  that  extensive 
experimental  work  with  porous  media  indicates  fluid  flow  is 
governed  by  Darcy's  law  for  the  range  of  Reynolds  numbers 
where  viscous  effects  dominate.   The  upper  limit  (velocity) 
Reynolds  number  in  these  investigations  is  subject  to  dis- 
agreement and  varies  from  0.1  to  75.   Inertial  effects 
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increase  with  increasing  Reynolds  numbers.   Boffa  [Ref.  12] 
has  shown  that  for  a  fixed  Reynolds  number,  inertial  effects 
diminish  with  increasing  air  temperature.   Darcy ' s  law  for 
two-dimensional  flow  is. 


a  =  -f.^^;.i^^.,,^,u)  ,3.5, 


where  Q  is  the  filter  velocity  or  the  volumetric  flow  rate 
per  unit  cross-sectional  area;  m  is  the  specific  permeability; 
g  is  the  gravitational  acceleration;  r  and  z  are  the 
unit  vectors  in  the  r  and  z  directions,  respectively;  and 
dP/dr  and  dP/dz  are  the  pressure  gradients  in  the  r  and  z 
directions,  respectively.   The  assumption  here  is  that  the  r 
(radial)  and  z  (axial)  velocity  components  each  react  to  the 
pressure  independent  of  the  other.   The  specific  permeability 
of  the  porous  medium  used  in  the  model  is. 


m   =   9^(7)^  (3.6) 


Expression  3.6  is  based  on  a  cappilaric-serial  model  given 
by  Scheidegger  [Ref.  10].  Expressions  for  permeability  vary 
with  the  physical  assumptions  of  the  flow  paths  incorporated 
into  the  model.  Permeability  also  varies  with  the  pore  size 
distribution  assumed.  Physically,  permeability  and  porosity 
are  not  related  [Ref.  9].  Porosity  is  a  quantifiable 
property  of  the  porous  medium,  whereas  permeability  is  a 
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constitutive  property  (i.e.,  a  property  specified  by  a 
constitutive  relation,  Darcy ' s  law).   The  specific  permea- 
bility is  proportional  to  the  filter  velocity  induced  by  a 
unit  pressure  gradient.   Darcy 's  law  is  not  derived  from 
first  principles  but  through  exhaustive  experimental 
analyses . 

The  Dupuit-Forcheimer  assumption,  addressed  in  Carman 
[Ref.  11],  relates  the  local  pore  velocity  to  the  filter 
velocity  by, 

Q   =   pV   ^  (3.7) 

The  hypothesis  is  that  the  local  pore  velocity  is  greater 
than  the  filter  velocity.   The  actual  velocity  in  a  single 
pore  is  a  function  of  the  fluid  element's  location  within 
the  pore.   The  Dupuit-Forcheimer  assumption  defines  an 
"average"  velocity  within  the  pore. 

The  continuity  equation  for  two-dimensional  flow  in 
porous  media  with  nonconstant  porosity  distributions  is, 


D(pp,) 

— =rTr=—  +  P  P  Div  V  =   0  (3.8) 

Dt     ^  "^a 


Invoking  the  Dupuit-Forcheimer  assumption  and  Darcy 's  law, 
the  continuity  equation  becomes. 


-3^—  +  pp^Div[— (3^  r  +  (3^  -  pp^g)z)]   =   0         (3.9) 
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Equation  3.9  (neglecting  body  forces)  is  one  of  four  field  equations 
cast  into  a  finite  element  formulation  later  in  this  work.   From  the 
pressure  distribution,  the  pore  velocity  distribution  is 
obtained  by  invoking  Darcy's  law  and  the  Dupuit-Forcheimer 
assumption.   The  derivation  of  Equation  3.9  is  presented  in 
Appendix  A. 

C.   SEMENOV  MODEL  OF  COMBUSTION 

The  model  of  Semenov  described  in  Frank-Kamenetskii 
[Ref.  13]  and  Vulis  [Ref.  14],  is  the  combustion  model 
employed  herein.   The  basis  of  the  model  is  the  relation  of 
reaction  rate  to  temperature  and  the  interaction  of  heat 
generation  and  heat  transfer.   The  reaction  rate  expression 
R  ,  is  the  Arrhenius  expression  for  a  simple  n-th  order 
reaction. 


R   =   A4)^exp( — —)  (3.10) 

^  R  T 

u  c 


where  A  is  the  time  constant  of  the  chemical  reaction,  E  is 
the  activation  energy,  R   is  the  Universal  Gas  Constant, 
T   is  the  absolute  carbon  temperature,  and  $  is  the  oxygen 
concentration.   In  a  simple  reaction,  the  reaction  rate 
depends  on  the  concentrations  of  reactants  and  not  on  the 
products.   The  heat  generated  by  the  reaction  is  obtained 
by  multiplying  Expression  3.10  by  the  heat  (enthalpy)  of 
combustion.   As  explained  in  [Ref.  13],  a  theory  of  combustion 
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can  be  constructed  only  if  certain  reasonable  assumptions 
are  made.   For  example,  in  order  to  regard  the  initial 
states  as  stationary,  one  must  neglect  the  reaction  rates 
at  these  states.   The  empirical  law  of  chemical  kinetics 
embodied  in  the  Arrhenius  expression  tells  us  that  the  reac- 
tion rate  never  goes  to  zero  but  falls  off  exponentially 
with  a  decrease  in  temperature.   Neglect  of  the  reaction 
rate  at  the  initial  states  is  necessary  to  achieve  initial 
equilibrium.   In  a  finite  range  of  temperatures  above  the 
initial  states,  neglect  of  reaction  rates  is  also  justifiable 

[Ref .  13] .   At  room  temperature  the  reaction  rates  are  on 

2 
the  order  of  l.E-15  Ibm-carbon/f t  -hr  or  smaller,  vice 

2 
l.E+3  Ibm-carbon/ft  -hr  at  1500  °F. 

D.   ARRHENIUS  LAW  OF  REACTION  RATES 

In  1934,  Parker  and  Hottel  [Ref.  15]  proposed  an 
Arrhenius  expression  for  the  reaction  rate  of  carbon  reacting 
in  air.   The  expression  assumed  a  simple  first  order  reaction 
for  the  combustion  of  carbon  in  air.   Frank-Kamenetskii 
[Ref.  13]  has  shown  that  the  Parker  and  Hottel  data  is  better 
correlated  by  a  fractional  order  reaction,  n  between  1/3 
and  2/3.   A  reaction  order  of  1/2  yields, 


R    =   2.065  X  10^  (R^  (}))"^^^exp(-57,240/R  -T  )  (3.11) 

2 

In  Expression  3.11,  R   is  in  units  of  Ibm-carbon/f t  -hr, 

R^      is  the  gas  constant  for  oxygen  (48.29  ft-lbf/lbm-R) , 
^2 

28 


(\)    is  in  Ibm/ft  ,  R   is  1.986  Btu/lbmole-R,  and  T   is  in 
Rankine.   In  order  to  determine  the  rate  of  heat  generation 
and  the  rate  of  oxygen  consumption,  the  chemical  reaction 
for  the  combustion  process  must  be  considered.   Although  many 
complex  chains  in  chemical  kinetics  generally  describe 
combustion,  a  simple  two-product  analysis  is  employed  in  this 
formulation.   For  nonvolatile  combustion  of  carbon  and 
oxygen,  two  reactions  that  describe  the  process  are, 


C   +   J  0^      -^      CO  (3.12) 


C  +  O2   ^   CO2  (3.13) 


where  0  denotes  oxygen;  and  C,  CO  and  CO^  denote  carbon, 
carbon  monoxide  and  carbon  dioxide,  respectively.   The  ratio 
of  the  mass  rates  of  carbon  monoxide  to  carbon  dioxide 
produced  increases  with  increasing  temperature.   Arthur 
[Ref.  16]  presents  an  expression  for  the  rate  ratio  as  a 
function  of  temperature  (in  Kelvin) . 


^°   =   2500  exp(-6240/T  )  (3.14) 


CO2  ^ c' 

The  expression  is  valid  for  temperatures  between  790  and 
1170  K  (1310-2110  degrees  Fahrenheit) .   As  a  result  of  this 
temperature  dependency,  the  stochiometric  ratio  and  the  heat 
of  reaction  are  functions  of  temperature.   Defining  the 
fraction  of  carbon  monoxide  being  produced  by / 
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^co  =  <^'/'^  ^  '^'>  '2-15) 


and  the  fraction  of  carbon  dioxide  as, 


^CO^   =   ^/'^  "  '^>'  '3.16 


the  heat  of  combustion  is  then  expressed  as, 


^"r   =   ^CO  ^"CO  ^  ^CO^  '«C02  <3-^^ 


Values  for  the  heats  of  combustion,  AH    and  AH    ,  as 
functions  of  temperature  may  be  obtained  from  JANAF  (Joint 
Army,  Navy,  and  Air  Force)  tables  [Ref.  17].   For  the  range 
of  temperatures  being  investigated  (80-2000  degrees  Fahren- 
heit) ,  the  heats  of  combustion  are  3966.3  Btu/lbm  for  F 
and  14,121.  Btu/lbm  for  F    .   Frank-Kamenetskii  [Ref.  13] 
points  out  that  over  narrow  ranges  of  temperature  it  is 
permissible  to  use  an  approximate  expression  which  correctly 
describes  the  reaction  rate.   The  stochiometric  ratio  (fuel 
to  oxygen)  of  the  overall  reaction  is. 


^R  =^cofco/'^co/co  ^  ^co^co^'        ''•"' 


where  f^„  is  the  stochiometric  ratio  for  the  reaction  3.12 
and  f  is  the  stochiometric  ratio  for  the  reaction  3.13. 
The  rate  of  heat  generation,  R  ,  and  the  rate  of  oxygen 
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consumption  may  now  be  expressed  by, 


Rg   =   AHjj  R_^  (3.19) 


Rq^   =   fR^  \  (3.20) 


Parker's  and  Hottel's  work  [Ref.  15]  and  Arthur's  work 
[Ref.  16]  were  conducted  with  specific  types  of  carbon. 
Tables  and  references  exist  in  Smoot  and  Pratt  [Ref.  18]  and 
Frank-Kamenetskii  [Ref.  13]  for  rate  expressions  using  other 
types  of  carbon.   The  present  model  allows  for  any  simple 
fractional  order  rate  expression  of  Arrhenius  type  to 
account  for  carbon  consumption  with  CO  and  CO-  byproducts. 
For  this  work  the  Parker  and  Hottel  rate  expression  as  modi- 
fied by  Frank-Kamenetskii  (n  =  1/2)  is  used. 

Particle  diameter  decreases  with  progressive  combustion. 
The  rate  of  decrease  depends  on  the  amount  of  carbon  con- 
sumed at  a  point  over  time.   Observations  have  shown  that 
the  effect  of  decreasing  particle  diameter  is  significant 
when  the  reaction  is  concentrated  in  a  small  region  of  the 
porous  medium.   To  account  for  this,  an  expression  for  the 
time  rate  of  diameter  change  as  a  function  of  reaction  rate 
is  derived  (Appendix  B) .   For  spherical  particles  the  equation 
is, 


d   =   -2R  Z  D^/(7Tp  d^)  (3.21) 
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where  p   is  the  bulk  mass  density  of  carbon.   The  diameter 
and  the  reaction  rate  are  functions  of  time  and  space. 

E.   CARBON  HEAT  TRANSFER  EQUATION  FOR  POROUS  MEDIA 

The  heat  transfer  equations  of  the  present  investigation  incor- 
porate:  (1)  radiation  (at  boundaries  where  applicable),  (2) internal 
convection,  (3)  conduction,  (4)  internal  combustion ,  (5)  temperature 
dependency  of  properties,  and   (6)  compressibility  effects 
of  air  into  a  two-dimensional  (cylindrical  coordinate) 
formulation.   The  carbon  energy  conservation  equation  is, 

8T 
(l-p)p^C^  -^   =   V-  (1-p)  (kg)  (VT^)  -hZ(T^-T^)  +  RgZ      (3.22) 

The  derivation  of  Equation  3.22  is  presented  in  Appendix  A. 
The  effective  conductivity,  k  ,  of  the  porous  solid  was 
proposed  by  Russel  [Ref.  24], 

p. 2/3  ,^,i.p,2/3, 

k    =   k^  — ^ (3.23) 

I  2/3     ,  ^   a  ,  -,     ,  2/3  ,   ,  , 
P      -P+]^(1-P      +P) 
c 

where  k   and  k   are  bulk  thermal  conductivities  of  carbon 
c      a 

and  air  respectively,  and  p'  =  1-p.  Russel 's  expression, 
which  is  based  on  an  electrical  analogy,  is  valid  for  the 
full  range  of  porosities,  0.0  to  1.0. 

Because  of  the  difficulties  encountered  in  a  radiative 
analogy  to  Fourier's  law  of  conduction,  the  particle  to 
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particle  radiative  exchange  in  the  porous  medium  is  omitted. 
The  difficulties  are:   (1)  the  geometry  does  not  easily  allow 
one  to  derive  an  expression  for  a  "sink"  temperature  to  use  in 
a  linearized  approximation  of  the  Stef an-Boltzman  equation,  and 
(2)  the  multi-v/avelength  characteristics  of  the  radiation 
phenomenon  are  not  easily  incorporated. 

F.   HEAT  TRANSFER  EQUATION  FOR  AIR  IN  POROUS  MEDIA 

The  internal  convective  heat  transfer  coefficient  of 
Yoshida,  Ramaswami ,  and  Hougen  [Ref.  33],  is  given  by, 


h   =   0.91  Re'  °'^-^[i|^C  G(C^ii/k^)  ^^^]f  (3.24) 

a    a    a       m 


where  ^    is  equal  to  1  for  spherical  particles,  G  is  a  pseudo 

mass  velocity  given  by  p  p  s  and  C   is  the  specific  heat 

a        a 

of  air  at  constant  pressure.   The  f   subscript  indicates 
properties  are  to  be  evaluated  at  film  temperature.   Re'  is 
a  pseudo  Reynolds  number  defined  by. 


Re'   =   G/z  p  t|;  (3.25) 

The  air  properties,  as  well  as  the  internal  convection  coeffi- 
cient, h,  are  temperature  dependent.   The  reaction  rate  in 
Equation  3.22  is  given  by  Expression  3.19. 

An  energy  balance  on  the  air  within  the  porous  medium 
obtains  the  second  heat  transfer  equation. 
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PPa^a  TF  =   V-<P^'T^)  *   hZ(T^-T^) 


-  pp  C^(V-V)T^  -  (V-V)pP  (3.26) 


The  density  of  air,  p  ,  is  approximated  by  the  ideal  gas  law. 

a 

Pressure  terms  are  due  to  the  compressibility  of  air.   The 
derivation  of  Equation  3.26  is  presented  in  Appendix  A.   All 
properties  in  Equation  3.26  are  temperature  dependent.   The 
properties  of  standard  air  were  used  in  the  model.   Vatikiotis 
[Ref.  9]  points  out  that  tolerable  differences  (average  of 
7%  difference)  between  standard  air  properties  and  properties 
accounting  for  the  presence  of  byproducts  CO  and  CO-,,  is 
acceptable.   Increased  accuracy  would  introduce  an  additional 
mass  balance  equation  for  either  CO  or  CO^.   Polynomial  ex- 
pressions used  to  calculate  the  thermophysical  properties  of 
air  are  those  presented  by  Vatikiotis  [Ref.  9].   The  expressions 
are  presented  in  Appendix  B. 

G.   OXYGEN  DIFFUSION  EQUATION  FOR  POROUS  MEDIA 

The  fourth  field  equation  necessary  to  complete  the 
system  of  equations  is  provided  by  an  oxygen  species  conser- 
vation requirement.   Transport  mechanisms  included  in  the 
model  are  molecular  diffusion  (Fick's  law),  convective  mass 
flow  and  oxygen  consumption  due  to  combustion.   Pressure  and 
temperature  induced  concentration  gradients  are  considered 
negligible.   Vatikiotis  [Ref.  9],  provides  examples  for  which 
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diffusion  arising  from  pressure  and  temperature  gradients 
is  important.   The  oxygen  mass  balance  equation  is, 


Pft   "   V-(pPgV4))-V-(p({)V)-RQZ  (3.27) 


The  derivation  of  Equation  3.27  is  presented  in  Appendix  A, 
The  effective  diffusivity,  V    ,    proposed  by  Denbigh  and 
Turner  [Ref.  31]  for  a  porous  medium  is 


V^      =      V/t  (3.28) 


Expression  3.28  accounts  for  the  tortuosity  encountered  by 
the  oxygen  molecules  as  they  flow  through  the  porous  medium. 
The  semi-empirical  expression  proopsed  by  Gilliland  [Ref.  34] 
is  used  to  obtain  the  diffusivity  of  oxygen  into  air.   The 
expression  is, 


^3/2   -1    -1  1/2      1/3    1/3  2 

V  =      435.7  T^/^(M^-'  +M_-')'-/^/[P(V^/-'  +V;!:/ ■')'']  (3.29) 

2 

where  V    is  in  units  of  cm  /s,  P  is  the  total  pressure  in  Pa, 

V  and  V    are  the  molecular  volumes  of  air  and  oxygen, 
a      O2 

respectively.   M   and  M    are  the  molecular  weights  of  air 

a      O-^ 

and  oxygen.   The  values  of  V   and  V   may  be  obtained  in 

a       (j-^ 

3 

Holman  [Ref.  35]  as  29.9  and  7.4  cm  ,  respectively.   The 

oxygen  consumption  term  in  Equation  3.27  is  given  by  Expression 
3.20. 
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H.   BOUNDARY  CONDITIONS 

Boundary  conditions  employed  were  as  follows  for  the 
carbon. 


3T 
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where  q   is  the  starting  heat  flux, 
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For  pressure. 
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Finally,  for  0_  concentration, 
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3.43a) 


P  V 


e     9z 


=  1 


v(4)  -4)^) 


3.44) 


P  V 


e     9r 


^  =  1 
r 
o 


(3.45: 


Equations  3.35,  3.36,  3.43,  and  3.44  correspond  to  convective 

flux  conditions  on  air  temperature  and  0^  concentration. 

At  the  air  inlet  ( —  =  0.0),  Dirichlet  conditions,  Equations 

o 
3.35a  and  3.43a  may  be  imposed  on  the  air  temperature  and 

oxygen  concentration.   The  idea  behind  this  treatment  is  that 

in  the  presence  of  a  semi-infinite  medium  (ambient  air) ,  the 

air  and  0^  concentration  may  be  considered,  to  a  first 

approximation,  very  near  ambient  conditions.   Equations  3.30, 

3.34,  3.38  and  3.4  2  correspond  to  symmetry  conditions  at 
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—  =  0.0.   Equations  3.33,  3.37,  3.41  and  3.45  correspond 

o 
to  impermeable  boundaries  and  insulated  boundaries  at  —  =  1.0 

o 
and  are  used  in  conjunction  with  one-dimensional  problems. 

The  second  part  of  Equations  3.33  and  3.37  are  the  boundary 

conditions  that  represent  a  relaxation  of  the  radial  insulation 

of  carbon  and  air  temperatures  which  allow  the  system  to 

become  a  two-dimensional  heat  transfer  problem.   It  must  be 

pointed  out  however  that  the  heat  transfer  coefficients  of 

Equations  3.33  and  3.37  are  not  easily  obtained.   This  work 

could  not  proceed  if  the  following  assumptions  were  not  made. 

It  was  assumed  that  the  air  and  carbon  had  near  equal  profiles 

at  —  =  1.0  so  that  the  same  heat  transfer  coefficient  would  apply 

o 
to  both  variables.   If  Equation  3.35  is  used  as  a  boundary 

condition,  the  air  temperature  follows  closely  the  carbon  tempera- 
ture at  —  =  1.0.   Furthermore  it  was  assumed  for  simplicity  of 

o 
calculation  that  the  heat  transfer  coefficient  was  constant. 

The  difficulty  is  as  follows.   In  order  to  determine  the 
heat  transfer  coefficient  for  the  impressed  temperature 

profiles  at  —  =  1.0,  the  temperatures  themselves'  must  be 

o 
known.   The  correct  procedure  would  involve  an  initial  esti- 
mate of  the  profile,  calculation  of  the  heat  transfer  coeffi- 
cient, solution  of  the  problem  for  one  integration  and 
verification  of  the  assumed  and  calculated  temperature 
profiles.   Upon  suitable  verification  of  the  profiles,  the 
same  procedure  is  performed  for  the  next  integration.   Other 
problems  arise  for  example  if  the  radial  heat  transfer  at 
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—  =  1.0  is  of  natural  convective  type.   At  the  exit  —  =  1.0, 

o  o 

functional  and  derivative  continuity  must  be  demonstrated 

for  pressure,  oxygen  concentration  and  air  temperature  at 
the  vertical  boundary  separating  the  stream  of  exiting  air 
and  the  convective  boundary  layer.   This  is  known  as  a  conju- 
gate problem.   The  above  method  is  a  realistic  method  for 
the  treatment  of  this  problem  but  admittedly  it  is  beyond  the 
scope  of  this  work.   The  numerical  code  has  provisions  for 
incorporating  an  average  of  isoflux  and  isothermal  natural 
convective  Nusselt  number  formulations  (Churchill  correlations 
for  vertical  cylinders  obtained  from  [Ref.  36].   Although  the 
subroutine  has  been  written  there  has  been  no  attempt  to 
date  to  implement  it.   It  is  known  that  physically  heat 
transfer  coefficients  lie  between  those  generated  from  iso- 
thermal and  isoflux  considerations.   However,  simply  having 
empirical  correlations  for  evaluating  the  heat  transfer 
coefficients  does  not  eliminate  the  need  of  the  iterative 
scheme  described  above.   Additional  considerations  must  be 
addressed  in  this  context.   In  order  to  effect  an  impermeable 

vertical  boundary  at  —  =  1.0,  the  system  must  be  enveloped 

o 
by  a  cylindrical  container  fabricated  of  some  type  of  material 

This  material  has  an  inherent  potential  to  absorb  internal 

energy  that  transits  from  the  initial  system  to  the  boundary 

where  some  flowing  medium  is  able  to  convect  energy  away 

(Figure  3.2).   The  material's  ability  to  absorb  energy  adds 

a  thermal  resistance  to  the  heat  transmission  in  an  electical 
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analog  sense  (Figure  3.3)  and  thus  affects  the  temperature 
profile  and  ultimately  the  heat  transfer  coefficients.   Of 
necessity  then,  the  material's  inherent  ability  to  absorb 
energy  is  denied.   No  attempt  is  made  here  to  conjecture  as 
to  how  physically  one  may  impose  the  radial  boundary  conditions 
implemented  below.   If  it  were  possible  to  rigorously  apply 
the  heat  transfer  boundary  condition  described  above,  the 
numerical  model  would  be  capable  of  obtaining  a  solution. 
The  model  was  tested  at  various  values  of  constant  heat  trans- 
fer coefficient:   0,  1,  5,  50,  100,  500,  to  determine  the 
effect,  of  the  heat  transfer  coefficient  on  the  solution.   It 
was  found  that  until  the  heat  transfer  coefficient  gets 
significantly  large  compared  to  the  start  flux  applied  to 

the  carbon  at  the  air  inlet  boundary  (in  this  case  above 

2 

100  Btu/ft  -hr  compared  to  1500  flux  applied  to  carbon) ,  the 

two-dimensional  effects  on  the  temperature  profile  were 

localized  in  the  —  =  0.85  to  1.0  region.   The  temperature 

o 
profiles  for  —  less  than  0.85  were  essentially  one-dimensional 

o 
(constant  value  independent  of  r) . 

I.   INITIAL  CONDITIONS 

The  model  is  developed  so  each  problem  begins  at  a 

uniform  initial  temperature.  A  heat  flux  applied  to  the 
carbon  along  a  boundary  provides  a  means  of  bringing  the 

system  to  a  temperature  level  where  reaction  terms  are 

sufficiently  high  to  generate  combustion.   The  heat  flux  is 
treated  as  follows. 


41 


(l-p)V(k  ) 


c 

e'  8z 


=   -q"  (3.46 

z 
o 


Expression  3.46  is  incorporated  into  the  numerical  formula- 
tion as  Neumann  boundary  condition  for  carbon  temperature. 
The  initial  heat  flux  may  be  turned  off  at  any  specified 
time.   The  boundary  condition  of  Expression  3.46  then 
switches  to  a  Cauchy  (mixed)  boundary  condition  to  account 
for  radiative  heat  transfer  from  the  carbon  particles  at  the 
boundary  to  the  ambient  air.   Convective  heat  transfer  from 
the  carbon  particles  to  the  air  is  treated  through  the 
internal  heat  transfer  coefficient.   The  above  procedures 
for  treating  problem  initiation  obviate  the  need  of  trying 
to  specify  initial  conditions  which  may  in  general  be 
arbitrary  for  each  problem. 
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Figure  3.1   Geometric  Model  of  Porous  Medium 
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Figure  3.2    Top  View  of  Radial  Convection  Circuit 
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Figure  3.3    Electrical  Analogy  of  Thermal  Circuit 
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IV.   NUMERICAL  CONSIDERATIONS 

A.   TREATMENT  OF  NONLINEAR  COMBUSTION  TERM 

The  treatment  of  the  nonlinear  combustion  term  is  now 
discussed  in  brief  detail  highlighting  the  advantages  of 
each  treatment.   A  more  comprehensive  discussion  is  presented 
in  Appendix  C.   The  interested  reader  is  encouraged  to 
review  the  detailed  analysis.   There  are  several  ways  to 
treat  the  Arrhenius  reaction  rate  expression. 
1.   Excitation  (Force)  Term  Treatment 

Perhaps  the  simplest  method  for  incorporating  the 
combustion  term  is  to  evaluate  it  at  each  time  step  and  use 
the  evaluated  value  (held  constant)  as  an  estimate  of  the 
mean  value  during  the  next  integration.   It  is  realized  that 
in  fact,  the  value  is  not  constant  in  real  time.   This  treat- 
ment, however,  provides  a  means  of  incorporating  the  term 
into  the  system  of  equations  as  a  first  approximation  with 
relatively  little  computational  effort.   Averaging  techniques 
exist  for  improving  upon  this  method.   In  this  scheme  the 
excitation  vector  is  modified  at  each  nodal  point  in  the 
carbon  temperature  and  oxygen  concentration  equations  to 
include  the  combustion  terms.   The  exponential  reaction  rate 
term  that  appears  in  both  the  porous  solid  and  oxygen  diffu- 
sion equation  is. 


R    =   A  (|)"exp(-E/RT  )  (4.1) 

C  T         JT  ^ 
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In  the  carbon  equation,  Expression  4.1  appears  as, 


Rg   =   AHj^  .  R^  (4.2 


In  the  oxygen  concentration  equation.  Expression  4.1  appears 
as , 


•^02   =   ^R^  •  ''c  <4-3' 


The  term  that  is  to  be  evaluated  at  the  last  time  step  and 
to  be  held  constant  over  the  next  time  step  is. 


*t 

R    =   {A  (|)^exp(-E/RT  )}   ^~  (4.4) 


where  the  superscript  *t  _,  indicates  evaluation  occurring 
at  the  previous  time  step. 

2 .   Linear  Operator  Treatment  of  0^  Concentration 

In  order  to  realize  an  improvement  over  the  first 
treatment  one  may  retain  a  portion  of  the  reaction  expression 
as  an  operator  by  making  the  following  rearrangement: 


*t 
R^   =   {A  4)^"^exp(-E/RT^)  }   ^""^  •  [(p]  (4.5 


where  [ip]    denotes  a  spatial  operator  treatment  of  the 
response  variable. 
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3.   Bilinear  Operator  Treatment  of  0^  and  T 
c 2 c 

The  bilinear  operator  treatment  of  the  reaction  rate 
is  the  present  method  of  treatment.   The  reaction  rate  term 
is  rearranged  as  follows, 


A  (f  exp(-E/RT  )  *t 

R    =   { —}      ^    ^    '     [T  ^]        (4.6 

^  T  ^ 

c 


Letting 


^c 


5^6^^  ,   i  =  1,3  (4.7) 


and 


^      =      5^64^  ,    i  =  1,3  (4.8 


and  invoking  natural  coordinates  [Ref.  28],  an  elementally 
averaged  contribution  results  in  the  following  area  integral, 

C   =   C  /  /   ?  c'^e,  c'^e   dA  (4.9) 

e 

The  double  subscript  permutation  of  carbon  temperature  and 
O2  concentration  results  in  a  3  x  9  elemental  matrix  that  is 
distributed  into  the  system  equations.   In  this  scheme  at 
each  nodal  point  within  an  element,  carbon  temperature  and 
O2  concentration  equations  receive  nine  terms  arising  from 
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combustion  occurring  within  the  element.   An  element  loop  is 
performed  to  account  for  combustion  terms  in  an  elementally 
averaged  sense  over  the  entire  system  spatial  domain. 

B.  OPTIMAL  COMPACT  STORAGE  (OCS) 

The  optimal  compact  storage  scheme,  presented  by  Franke 
and  Salinas  [Ref.  29],  employed  in  the  solution  of  the 
combustion  problem,  allowed  the  computational  effort  to 
proceed  with  substantial  savings  in  computer  storage, 
computer-run  time,  and  ultimately  dollar  cost  per  run.   In 
the  optimal  compact  storage  scheme,  only  the  non-zero  coeffi- 
cients are  stored  in  a  vector  array.   This  results  in  a 
very  significant  reduction  of  the  storage  area  compared  to 
banded  storage.   One  might  encounter  problems  on  the  order 
of  1000  X 1000  DOF.   (In  this  problem  there  are  four  degrees 
of  freedom  at  each  nodal  point.   A  17  x 17  grid  requires  a 
1156  X 1156  matrix  if  full  storage  is  employed.)   For  banded 
storage  the  bandwidth  might  be  200.   The  storage  ratio  for 
bandwidth  to  full  storage  would  be  0.2.   Using  OCS,  a  con- 
servative estimate  of  the  average  number  of  equation  entries 
in  this  model  is  20.   The  ratio  of  OCS  to  banded  storage 
is  0.1.   Thus  in  this  example,  the  savings  realized  by  OCS 
vice  full  storage  is  98%! 

C.  GRID  CONVERGENCE 

A  convergence  study  was  done  on  the  response  field  for 
several  grids.   It  was  found  that  non-uniform  grids  were  far 
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superior  to  uniform  grids.   In  fact,  for  many  problems, 
computer  storage  limitations  require  that  a  non-uniform  grid 
be  used  to  obtain  reasonably  accurate  solutions.   Non- 
uniform grids  allow  the  investigator  to  take  advantage  of 
the  knowledge  of  areas  of  severe  combustion  induced  gradients. 
"Stacking"  elements  in  these  areas  assists  the  algorithm  in 
producing  more  accurate  and  numerically  stable  results.   For 
similar  degrees  of  accuracy  with  uniform  grids,  it  is  esti- 
mated that  grids  on  the  order  of  1000  x 1000  nodal  points 
(and  larger)  would  have  been  required.   Nondimensionalized 
length  discretizations  on  the  order  of  .001-. 005  and  smaller 
near  the  air  inlet  boundary  permit  excitations  on  the  system 
boundaries  to  be  propagated  accurately  through  the  medium. 
Until  discretizations  on  the  order  of  non-dimensionalized 
lengths  of  0.001  were  employed,  numerical  instability  was 
encountered  in  the  carbon  temperature  field.   This  instability 
was  manifested  by  severe  overshoot,  leading  to  negative 
temperatures  in  close  proximity  to  severe  thermal  gradients. 

For  non-uniform  grids  it  was  found  that  grids  below  10 
nodal  points  in  each  direction  were  not  adequate.   The 
temperatures  obtained  from  these  models  were  lower  due  to  the 
relatively  low  degree  of  grid  refinement  or  discretization 
possible  with  such  few  points.   Grid  refinement  in  this  type 
of  problem  is  essential  to  capture  the  high  activity  (large 
gradients)  that  is  inherent  in  the  combustion  phenomenon. 

A  convergence  study  was  done  on  several  nodal  points  at 
various  times  for  non-uniform  grids  at  12  xl2,  17  xi7,  and 
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22  X 22  nodal  points.   The  solution  changed  5%  from  a  12  x 12 
to  a  22  X 22  nodal  point  model.   Next,  a  17  x 17  point  non- 
uniform grid  model  was  investigated.   It  was  found  that 
the  solution  changed  approximately  2%  from  the  22  x  22  nodal 
point  model.   At  this  point,  it  was  determined  that  the 
12  X 12  non-uniform  nodal  point  afforded  the  desired  balance 
between  cost,  computational  effort  and  computer-run  time. 


grid  size 

kmax 

cloc 

k 

runtime 

(min. ) 

(n  X 

n) 

12 

6,398 

70 

17 

13,223 

100 

22 

22,498 

110 

D.   MODEL  VALIDATION  TESTS 

In  order  to  validate  the  capabilities  and  accuracy  of 
the  numerical  code,  several  tests  were  conducted. 

1.   The  Steady  State  Problem 

First,  it  was  necessary  to  ensure  the  model's  ability 
to  recognize  a  steady  state  condition.   This  was  the  simplest 
of  all  tests  (and  first  to  be)  performed.   Initial  (ambient) 
conditions  were  imposed  uniformly  throughout  the  system  on 
the  four  fields:   carbon  temperature  and  air  temperature 

(80  degrees  F.),  pressure  (2,116.8  psf)  and  oxygen  concen- 

3 
tration  (0.0172  Ibm/ft  ).   An  initial  integration  time  step 

equal  to  the  minimum  time  step  allowable  (user  input)  l.E-6  hr . , 

was  selected.   In  less  than  five  integrations  the  algorithm 

switched  to  the  maximum  allowable  time  step  (user  input) 
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5.E-1  hr.,  and  iterated  about  the  given  initial  fields. 
This  sufficiently  demonstrated  the  first  of  several  tests 
requisite  to  assure  numerical  code  compatibility  with  the 
integration  scheme. 

2.   The  Finite  Slab 

Consider  the  slab  of  thickness  2L  in  Figure  4.1(a) 
One  seeks  the  solution  of  the  one-dimensional  conduction 
equation , 


T—     =      a   — ^  (4.10) 

8x 


subject  to  an  initial  condition, 


T(x,0)   =   T^  (4.11) 


and  uniform  convection  conditions  at  both  surfaces:   At 
X  =  ±L: 


-k  1^   =   ±  h(T  -T  )  (4.12 

dX  O 


The  exact  solution  is  given  as  a  Fourier  series  in  [Ref.  23], 


l^     =        I      C.  e-e'^t/L^cosCS.  |)  (4.13) 

i    o      i=l 


where. 


52 


4  sin  6  • 
"^i   "   23.  +sin(26.)  (4.14, 


The  constants  S ■  are  the  roots  of  the  transcendental  alge- 


braic equation, 


!.  tan  3.   =   B.   =   h  L/k  (4.15 

1       1        1        o 


Thus  the  solutions  have  the  Biot  number  of  the  slab  as  a 
parameter.   The  roots  of  Equation  4.15  are  tabulated  in 
Appendix  A  to  [Ref.  25]. 
3.   The  Cylinder 

Consider  the  cylinder  of  Figure  4.1(b).   The  suddenly 
immersed  cylinder  satisfies  the  equation^ 


3T      a  3  ,   3T,  ,,  ,,. 

It  =  F  37^^  37^  (^-^^^ 


subject  to  an  initial  condition, 


T(r,0)   =   T^  (4.17) 


and  a  convection  condition  at  the  surface:   At  r  =  r  , 

o 


-k  1^   =   h  (T  -T^)  (4.18) 

3r      o     o 


The  solution  is  given  as  a  Fourier  series  in  [Ref.  23] : 
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2     2 

T  -  T         °°      -3  •  at/r 

1     O        1=1 


where. 


2      ^l^^i^ 
C    =   ^  -^ :^ ^ (4.20 


The  constants  3'  are  the  roots  of  the  algebriac  equation, 


B.J^(B.)/J^(6.)   =   B.   =   h^r^k  (4.21) 


The  functions  J   and  J,  are  the  Bessel  functions  of  the  first 

o       1 

kind  whose  numerical  values  are  tabulated  in  most  advanced 

engineering  mathematics  texts.   For  a  limited  range  of 

interest,  numerical  values  of  J   and  J,  are  tabulated  in 
'  o       1 

Appendix  A  of  [Ref.  25].   The  roots  of  Equation  4.21  are 
tabulated  in  [Ref.  23]. 

4 .   Multidimensional  Solutions  by  the  Product  Method 

The  classic  solutions  for  semi-infinite-  and  finite- 
thickness  bodies  (slabs,  cylinders  and  spheres)  may  be  used 
to  generate  solutions  for  multidimensional  bodies. 

Use  of  the  product  solution  technique  is  illustrated 
by  Figure  4.2.   The  problem  chosen  to  validate  the  heat 
transfer  portion  of  the  current  model  is  the  "sudden  immer- 
sion" problem  for  a  finite-length  cylinder  subjected  to 
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uniform  convection  conditions  (h  ,T  )  on  all  sides  and 

o  °° 

initial  uniform  temperature  T..   The  differential  equation 
to  be  solved  is. 


19,  ^e^  ^  3^e  i  ae             .,  _^ 

-  TT—  r  T—   +  y  =   -  TT                 4.22 

r  8r    9r     ^  2  a  3t 

dX 


where  9  =  (T-T  )/(T.  -T  ).   The  initial  condition  is 


e(r,x,0)   =   1  (4.23) 


The  convection  boundary  conditions  are, 


At  top  and  bottom:      -k  4^  =   ±  h  6  (4.24) 

^  8x         o 


On  the  sides:  -k  |^  =   he  (4.25) 

dr      o 


The  solution  is  the  product  of  the  two  simpler  analyses: 
the  semi-infinite  slab  and  the  infinite  cylinder.   Let 


(x,r,t)   =   ®slab^^'^^  *  ^cyl^^'"^^  (4.26) 


=   P(x,t)  •  C(r,t) 

Substituting  Equation  4.26  into  Equations  4.22,  4.23  and 

4.24  and  separating  the  variables,  reduces  the  two-dimensional 

problem  into  two  one-dimensional  problems: 
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Slab 


a^e 


3x 


1   s 

a   9t 


(6   =  0  -  ,  ) 
s     slab 


(x,0)   =   1 


30 


-k 


8x 


X  =  ±L 


h^03(±L,t) 


1  3     ^^c 
cylinder:   -  ^(r  -^; 


,  30 
1  c 

a   3t 


(0   =0   ,  ) 
c     cyl 


(r,0)   =   1 


30 


-k 


3r 


h  0  (r  ,t) 
o  c   o 


Thus  Equation  4.26  is  the  exact  solution  to  the  sudden 
immersion  problem  of  a  right  circular  cylinder. 
5 .   Validation  Problem 

The  validation  problem  is  Example  4.9  in  [Ref.  25]. 
The  problem  is  stated  here  for  the  convenience  of  the  reader. 

The  short  cylinder  in  Figure  4.3  is  initially  at  40°  C 

and  then  plunged  into  a  fluid  with  h  =  300  W/m-K  and  T  =  200°  C 

-9   2 
The  material  is  bronze,  k  =  26  W/m-K  and  a  =  8.6  •  10    m  /s. 

Find  the  temperature  at  the  center  of  the  cylinder  after 

5  minutes. 
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Solution:   176°  C  (348.8°  F) .   The  details  of  the  solution 
are  found  on  pages  186-187  of  [Ref.  25].   The  product  solution 

technique  was  utilized  in  obtaining  a  5-term  approximation 

2 
of  the  series  solution  for  Fourier  moduli  (F   =  at/L  , 

equivalent  to  a  nondimensionalized  time)  less  than  0.2.   For 

Fourier  moduli  greater  than  0.2,  a  two-term  approximation  of 

the  exact  solution  was  used.   Heisler  [Ref.  26]  points  out 

that  a  one-term  approximation  of  the  exact  series  solution 

yields  an  accuracy  to  within  1%  of  the  exact  value  for 

Fourier  moduli  greater  than  0.2.   Various  values  (5,  10,  20 

and  30  seconds)  corresponding  to  values  of  Fourier  moduli 

less  than  0.2  were  used  in  comparing  model  solutions  with 

the  exact  (five -term  approximations  to  the  exact)  solution. 

Values  (3,  5,  7  and  9  minutes)  corresponding  to  Fourier 

moduli  greater  than  0.2  were  used  to  compare  model  solutions 

with  the  exact  (one-tenn  approximations  with  1%  accuracy) 

solutions.   Figures  4.4  through  4.8  show:   (1)  the  effects 

of  time  on  the  solution  of  various  grids,  and  (2)  effects 

of  grid  refinement  at  equal  values  of  time.   Values  plotted 

on  the  ordinate  scales  are  deviation  (percentage  error)  from 

the  exact  solutions  versus  time.   Observation  yields  that  for 

small  values  of  time,  grid  refinement  obtains  more  accurate 

results  and  obtains  a  faster  convergence  to  steady  state. 

6.   Testing  the  Model's  Ability  to  Accept  and  Reject 
Heat 

An  equivalent  simple  test  in  model  validation  is  the 

model's  ability  to  accept  heat  from  a  boundary  flux  and 
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reject  it  at  a  later  time  to  its  environment  via  convection. 

The  applied  heat  flux  at  —  =  0  is  1500  Btu/ft  -hr. 

o 
Figures  4.9-4.12  illustrate  system  response.   Figure  4.13 

is  the  input  data  set  used. 

7 .   The  One-Dimensional  Problem 

Another  step  in  the  validation  of  the  two-dimensional 
problem  is  an  examination  of  a  one-dimensional  problem.   One 
might  infer  that  a  two-dimensional  model  includes  the  one- 
dimensional  model  as  a  subset.   This  was  demonstrated  in  a 
separate  test.   A  one-dimensional  problem  may  be  imposed  on 
the  model  in  two  ways.   One  way  is  to  make  the  length  to 
diameter  ratio  very  large.   (No  restrictions  on  excitations 
is  implied. ) 

Another  method  for  eliciting  one-dimensional  behavior 
is  to  "excite  the  system  in  a  one-dimensional  fashion," 

2  Z 

i.e.,  apply  boundary  conditions  at  —  =  0.0  and  —  =  1.0 

o  o 

independent  of  radius  and  insulate  radial  boundaries  at 

—  =  0.0  and  —  =  1.0.   The  insulated  boundary  at  —  =  0.0 
r  r  -^     r 

o  o  o 

arises  from  an  axisymmetric  formulation  of  the  problem.   The 

zero  gradient  (completely  insulated)  boundary  at  —  =  1.0 

o 
corresponds  to  a  completely  (all  four  variables)  impermeable 

vertical  boundary.   For  this  problem,  an  adiabatic  restriction 

is  sufficient.   In  general,  an  adiabatic  condition  means  an 

insulation  of  heat.   In  this  problem,  the  implication  is  far 

greater  because  of  the  coupling  between  system  fields.   An 

adiabatic  restriction  at  —  =  1.0  prohibits  non-zero  radial 

^o 
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pressure  gradients.   Unconstrained  radial  pressure  gradients 

at  —  =  1.0  imply  radial  convection  of  air  enthalpy  and  in  a 

o 
non-isothermal  fluid,  this  convection  undermines  the  initial 

adiabatic  assumption  at  —  =  1.0.   Oxygen  gradients  are  also 

o 
negated  by  an  adiabatic  condition  since  oxygen  convection 

by  air  must  occur  under  the  presence  of  nonzero  pressure 

gradients  at  —  =  1.0. 
^  r 

o 

The  two  methods  described  above  distinguish  one-dimen- 
sional behavior  arising  from  geometrical  conditions  and  one- 
dimensional  behavior  resulting  from  restrictions  on  excitations 
Figures  4.14-4.16  depict  model  one-dimensional  solutions  to 
"one  dimensional  excitations."   Figure  4.17  is  the  input  data 
set  used. 
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Figure  4.1    Classical  Transient  1-D  Geometries 
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Product  solution: 

eix.  r.  t)  =  P{x,  t)C(r,  t) 


Figure  4.2    Illustration  of  the  Product  Technique 
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Figure  4.3    Sample  Problem,  Cylinder  Geometry 
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Figure  4.4    Grid  Effects  at  Point  A,  (0,0: 
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Figure  4.5    Grid  Effects  at  Point  B,  (0,1) 
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Figure  4.6    Grid  Effects  at  Point  C,  (1,1) 


65 


TIME  (SECONDS) 


Figure  4.7    Grid  Effects  at  Point  D,  (1,0 
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Figure  4.8    Grid  Effects  at  Point  E,  (.5,.  5; 
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Model   Accepting   Heat    (Carbon 
Temperature   Profile) 
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Figure   4.10      Model   Rejecting   Heat    (Carbon 
Temperature   Profile) 
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Figure    4.11      Model    Accepting   Heat    (Oxygen 
Concentration   Profile) 
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Figure    4.12 


Model   Rejecting   Heat    (Oxygen 
Concentration   Profile) 
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Figure    4.13        Accept/Reject   Heat   Problem   Input   Data   Set 
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one  Dimensional  Problem  (Carbon 
Temperature  Profile  and  Oxygen 
Concentration  Profile) 
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Figure    4.15 


One   Dimensional   Problem    (Carbon 
Temperature   Profile) 


74 


03COHCt«TR»no<«suiir*a 

PflOaUUTlMliS    «Hi-w    M« 
„Otnu«»W"0-    'fOO    3iU-.T2HH 

no  -  iQ''o'  fT 
zo  -  lo'io'  rr 
OJ  conc  mm  -    "itt 


0!  c(»aNT««Tio  sb«f»a 

«0«LLM  TIME  IS    9C»-10'M>»t 

Kt»rf\u«»T  zyzo-  1600  snj/rrj-M^ 

no-    lO-lo'  fT 

ZO  -    1010*  FT 
03  CONC  MiH  -    0 


01  co»<iwTn*io«  SUM  *a 

woaUM  n-*  a  ui-o'  "*     ^  ^ 

no  -   lono'  FT. 
20  -    1  OIO'  FT. 

ojcoNCMM-  a 


03  CWCO<nWT10«  SUKF4CI 

moeuvi  nut  a  !X!-«'  -"* 
^tftux.Tivzo-  'MO  im/rn--* 

no-  lOno"  FT 
20  -  lO^lo'  FT 
03  COK  "IN  -    0 


Figure   4.16 


One   Dimensional   Problem    (Oxygen 
Concentration   Profile) 
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V.   DISCUSSION  AND  CONCLUSIONS 

Field  problem  behavior  depends  on  three  major  factors: 
(1)  Boundary  conditions,  (2)  System  excitation,  and  (3)  The 
physical  characteristics  of  the  system. 

Time  limitations  did  not  permit  this  investigator  to 
conduct  an  exhaustive  analysis  of  each  of  these  factors. 
Some  preliminary  analyses  were  performed  to  obtain  some 
understanding  of  the  effect  of  boundary  conditions,  and 
excitation  on  system  behavior. 

There  are  a  variety  of  physical  parameters  which  govern 
system  response,  such  as  permeability,  porosity,  and  the 
cylinder  length-to-diameter  ratio.   Here  only  a  brief 
investigation  of  the  effect  of  porosity  on  system  behavior 
was  undertaken,  and  is  reported. 

A.   BOUNDARY  CONDITIONS 

The  boundary  conditions  used  in  the  present  investiga- 
tion were  presented  in  Chapter  III.   Other  boundary  conditions 
are  possible.   For  example,  if  in  place  of  Ecp .    3.32  and  3.35,  an 

insulated  boundary  condition  at  —  =  1.0  is  imposed  on  the 

o 
carbon  and  air  temperatures,  then  preliminary  results  indi- 
cate temperature  response  is  higher  for  equal  values  of  time. 
This  behavior  is  due  to  the  buildup  and  storage  of  energy 
within  the  system  associated  with  an  insulated  boundary. 
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Similarly,  in  place  of  Equation  3.44,  a  zero  flux  boundary 

condition  at  —  =  1.0  could  have  been  imposed  on  the  oxygen 

o 
concentration.   The  zero  flux  condition  implies  an  impermeable 

boundary  with  respect  to  oxygen  concentration  fluxes.   This 

boundary  condition  would  apply  to  a  very  long  cylinder  (i.e., 

a  regularity  condition).   In  this  investigation,  a  Cauchy 

(convective  flux)  boundary  condition  on  the  oxygen  concentra- 

tion  is  imposed  at  —  =  0.0  and  1.0.   As  oxygen  is  locally 

o 
consumed  in  the  interior  of  the  system,  oxygen  gradients  are 

created.   According  to  Pick's  Law  of  Diffusion,  a  depleted 

oxygen  region  is  replenished  by  the  diffusion  of  oxygen  from 

regions  of  relatively  high  concentration  to  regions  of  low 

concentration.   Thus  a  convective  flux  boundary  condition 

on  oxygen  causes  the  depletion  of  oxygen  concentration  to  be 

retarded  through  the  diffusion  mechanism. 

B.  EXCITATIONS 

The  effects  of  high  and  low  values  of  the  heat  flux 

2 

applied  at  —  =  0.0,  are  predictable.   High  values  of  heat 

o 
flux  lead  to  an  accelerated  development  of  both  the  carbon 

and  air  temperature  profiles  and  0^  concentration  depletion 

2 

within  the  system.   For  high  values  of  heat  flux  at  —  =  0.0, 

^o 

the  numerical  integration  scheme  eventually  slows  noticeably 
as  a  result  of  steep  gradients  observed  at  this  boundary. 

C.  PHYSICAL  CHARACTERISTICS  OF  THE  SYSTEM 

The  physical  characteristics  of  the  system  generate  a 
significant  effect  in  field  problem  solutions.   In  the 
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problems  discussed  here  the  transport  properties  k  ,  k  , 

e    3, 

'm/u,    V      and  the  Arrhenius  reaction  rate  term  are  governing 
factors  in  the  respective  field  equations  in  which  they 
appear.   As  previously  discussed,  the  net  balance  between 
heat  transfer  rates  and  heat  generation  dictate  whether  a 
reaction  will  proceed  to  extinguishment  or  combustion. 
There  are  two  time  constants  to  be  considered  in  a  combustion 
problem  of  this  type:   a  time  constant  for  the  Arrhenius 
reaction  and  a  time  constant  associated  with  momentum 
transport.   The  momentum  time  constant  affects  the  rate  at 
which  heat  is  convected  (removed)  out  of  a  differential 
volume.   The  reaction  time  constant  affects  the  rate  at 
which  heat  is  generated  (added  into)  within  a  differential 
volume . 

D.   POROSITY  ANALYSIS 

In  this  problem  there  are  many  parameters  one  might  vary 
in  order  to  examine  resultant  system  behavior.   Due  to  time 
limitations  this  investigation  examines  the  effect  of  one 
parameter,  porosity,  on  system  behavior.   The  observed 
effect  of  porosity  values  ranging  from  0.476  to  0.90  is 
discussed.   Of  the  nine  runs  attempted,  only  five  had  suffi- 
cient output  data  that  allowed  for  comparative  analyses. 
Employing  Equation  3.1  for  the  porosity  associated  with 
spherical  particles,  the  carbon  diameter,  d,  was  varied  to 
achieve  various  values  of  porosity  while  holding  the  parti- 
cle center-to-center  distance  D  (Figure  3.1),  fixed  at 
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0.000417  ft.  (0.005  in.).   Values  of  d  (ft.)  equal  to 
0.000417,  0.000411,  0.000396,  0.000385  and  0.000370  were 
used  to  give  values  of  porosity  of  0.476,  0.5,  0.55,  0.60, 
and  0.65,  respectively.   The  results  are  graphically  pre- 
sented in  Figures  5.1--5.4.   Figures  5.1  and  5.2  show  the 
carbon  temperature  profile  and  the  oxygen  concentration 
profiles  for  porosity  p,  =  0.476  (D  =  0.000417  ft.  =  0.005 
in.).   In  all  cases,  excitations  for  these  runs  are  1500 

Btu/ft  -hr  (in)  at  —  =  0 . 0  and  50  Btu/ft^-hr  (out)  at 

o 

r  z 

—  =  1.0;  the  pressure  is  14.65  psi  at  —  =  0.0  and  14.55 
r  ^  ^       z 

o  o 

at  —  =  1.0;  the  cylinder  length  to  diameter  ratio  is  0.5 

o 
(length  =  1.0  ft.).   Figures  5.1  and  5.2  represent  typical 

graphical  results  for  carbon  temperature  and  oxygen  concen- 
tration profiles.   Figures  5.3  and  5.4  are  a  tabular  summary 
of  the  results  of  the  carbon  temperature  profile  and  the 
oxygen  concentration  profiles,  respectively,  for  varying 
porosities.   Profiles  with  higher  values  of  initial  porosity 
are  observed  to  have  accelerated  development  of  temperature 
profiles  (i.e.,  as  the  fuel  diameters  decrease  (increasing 
porosity) ,  the  carbon  temperature  responds  at  a  faster  rate 
to  system  excitations) .   The  oxygen  concentration  response 
is  as  expected  (i.e.,  as  reaction  rate  goes  up,  the  oxygen 
concentration  goes  down) .   In  the  runs  made  for  porosity 
values  of  0.70,  0.75,  0.85  and  0.90,  numerical  difficulties 
were  encountered.   The  temperature  profiles  were  observed  to 
increase  sharply  compared  to  lower  values  of  porosity  (300-400 
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2 

F  in  the  first  minute)  at  the  —  =  0.0  boundary.   Due  to 

z  -^ 

o 

the  steep  gradients,  further  system  behavior  could  not  be 
analyzed  since  output  data  consisted  of  a  single  output 
(T  =  5-20  S.,  problem  time,  depending  on  porosity  value) 
in  15  minutes  of  CPU  time. 

E.  CONCLUSIONS 

The  combustion  analysis  program  (CAP)  is  a  viable  tool 
for  the  analysis  of  heat  transfer  in  porous  media.   The 
model  constructed  provides  the  user  with  the  flexibility  to 
solve  problems  with  or  without  combustion.   The  role  of 
the  permeability  model  must  not  be  underestimated.   For  it 
is  the  physical  assumptions  in  this  aspect  of  the  model  that 
govern  the  flow,  heat  transfer  and  in  the  end  the  evolution 
of  the  combustion  problem. 

F.  RECOMMENDATIONS 

Follow-on  work  is  recommended  in  the  following  areas: 

1.  Develop  a  restart  capability  for  the  program  which 
will  allow  for  a  problem  to  begin  at  any  point  in 
time  with  initial  conditions  and  rate  information 
being  identical  to  previous  timestep  values. 

2.  Nondimensionalize  equations  in  terms  of  other  system 
quantities  in  addition  to  spatial  dimensions. 

3.  Flex  the  model  under  other  system  parameters  to  examine 
effects  on  system  behavior. 


4.  Consider  an  internal  (particle-to-particle)  radiative 
heat  transfer  analysis. 

5.  Consider  other  fuels  and  more  detailed  chemical 
kinetics  chains. 
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Figure    5.1 
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Figure    5.2 


Oxygen   Concentration 
Radial    Heat   Transfer 


(Porosity    =    0.476) 
(Out)     50    Btu/ft^-hr 
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Figure  5.3    Maximum  Carbon  Temperature  Summary  Porosity: 
p^  =  0.476,  p2  =  0.5,  p   =  0.55,  p   =  0.60,' 
P5  =  0.65.   Time:   t,  =0.1  minute? 
t^  =  1  minute,  t^  =  2  minutes,  t.  =  4  minutes 
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Figure  5.4    Oxygen  Concentration  Profile  Summary 

Porosity:   p,  =  .476,  P2  =  0.5,  p..  =  0.55, 
p^  =  0.60,  p^  =  0.65.   Time:   t,  =  0.1 
minute,  t-  =  1  minute,  t-.  =  2  minutes, 
t ,  =  4  minutes 
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APPENDIX  A 
FORMULATION  OF  FIELD  EQUATIONS 

A.   PRESSURE  DISTRIBUTION  EQUATION 

Darcy's  law  for  two-dimensional  flow  is, 


where  Q   and  Q   are  expressed  as  follows, 


Q    =   -  -  I?  (A. 2) 

r        p  9r 


^         m , 8P         ,  / ,  T 


Invoking  the  Dupuit-Forcheimer  relation,  and  solving  for  the 
pore  velocity  components  u  (radial  velocity)  and  v  (axial 
velocity).  Equation  A.l  becomes. 


u   =   —  ^—  A. 4 


~ni ,  8P         s  /  -^  r  \ 


The  continuity  relation  (derived  in  Appendix  B)  is. 


D(pp_) 

— ^^^-  +  p  p^Div  V   =   0  (A. 6) 
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Substituting  Equation  A. 2  and  Equation  A, 3  into  Equation  A. 6 
yields, 


^  F^'H^  ^  'H-  PP^g)^^)-Vpp,  (A. 7) 


Expanding  terms,  Equation  A. 7  becomes. 


.,    2        7    1  p      3r  m8r        y8r        r8r 

dr  dz  a 


^     ,  1     ^^a   ^    1    9m         1    3y,   .9?    ^  .  y       ^Pp 

p,    dz  m    9z         p    9z'  '9z         ^  ^a^         p    m      9t 

a  "^a 


a 


(A. 


Equation  A. 7  with  associated  boundary  conditions  (presented 
in  Section  II. B)  is  cast  into  a  finite  element  formulation 
and  becomes  one  of  four  field  equations.   Pressure  gradient 
information  is  then  substituted  into  Equation  A. 4  and  Equation 

A.  5  to  obtain  the  pore  velocity  components. 

B.  POROUS  SOLID  HEAT  TRANSFER  EQUATION 

In  performing  energy  balances  on  both  the  porous  solid 

and  on  the  air,  a  differential  volume  of  porous  medium  may 

be  partitioned  into  respective  volumes,  dV   =  (l-p)dV  for 

the  solid,  and  dV   =  pdV  for  the  air  (shown  in  Figure  A.l) . 

a 

The  convention  used  for  the  energy  balance  of  an  arbitrary 
differential  volume,  dV,  is. 
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Heat  into  DV  +  Heat  generation   =   Heat  out  of  DV 

+  Increase  in  internal  energy     (A. 9) 

The  heat  transfer  mechanisms  considered  for  the  carbon 
are  conduction,  radiation  heat  transfer  between  particles, 
convection  heat  transfer  from  the  particles  to  the  air, 
and  heat  generation.   Applying  the  above  convention,  the 
energy  balance  on  a  differential  volume  of  porous  solid  is 
(invoking  Tayler  Series  expansions  and  neglecting  higher 
order  terms) , 


[-  ?|F((l-P)rq^ond,r'  "  ll' '1"?' <3cond,  z' 


i(|7((l-p)rq^^d,^)  -  |l((l-p)q^^^_J]dV 


z 


-  q     dA'  +  q    dA'   =   (l-p)q.  ^dV 
^conv      ^gen  '^  ^int 


In  vector  form.  Equation  A. 10  becomes. 


-V-(q    ^  +  q   ^)dV  -  q     dA'  +  q    dA'   =   (l-p)q.  ^dV 
^cond    ^rad      ^conv       gen  int 


(A. 10) 


(A. 11) 


Substituting  the  following  expressions  into  Equation  A. 10, 


q    -,   =   -  k  VT  Fourier's  law  (A.  12) 

^cond        e   c 
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q   ,    =   -  k  VT 
^rad         r   c 


Fourier  analogy 
(radiative) 


(A. 13 


q       =   h(T  -T  )       Newton's  Cooling  Law      (A.  14 
conv        c    a  ^j  v 


'gen 


=   R 


^int    =   Pc  ^c  I? 


Heat  generation 


Internal  energy 


;a.15) 


(A. 16: 


yields , 


lr9 


9T 


3T 


±(|_(,(l_p)(k^)  _^)  ,  |_((i_p)(k^)  _^)Hv 


-h(T  -T)dA'+RdA'   =   (l-p)p  C  ^  dV 
c    a        g  '^   c  c  3t 


:a.i7 


Dividing  through  by  dV,  and  defining  dA ' /dV  as  Z,  the  specific 
internal  area  (i.e.,  surface  area  per  unit  volume),  Equation 
A.  17  becomes. 


9T 


V-  (  (1-p)  (k^)  VT^)  -hZ(T  -T)  +  R^Z   =   (l-p)p^C^^   (A. 18) 
ec        ca     g  ccdt 


The  expressions  used  to  obtain  values  of  the  properties  and 
parameters  in  Equation  A. 18  are  presented  in  Section  III.E. 

C.   AIR  HEAT  TRANSFER  EQUATION 

The  formulation  of  the  air  heat  transfer  equation  begins 
with  the  general  two-dimensional  energy  balance  equation, 


The  difficulty  in  obtaining  an  expression  for  k   is 
addressed  in  Section  III.E.   Throughout  this  work,  one  may 
keep  in  mind  that  kg  should  really  be  (kg  +kj-)  to  account 
for  conduction  as  well  as  radiation  within  the  porous  medium, 
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p  p  D(U  +K) 
a 

Dt 


->-->-     ->■->- 


=   -V-pq  -  V-pPV  -  (V-(pT-V) 


+  hZ(T^  -T^)  +  p  p^  (V-g) 

Co.  3. 


(A. 19) 


where  U  is  the  internal  energy,  K  is  kinetic  energy,  t  is 

the  dissipation  function,  and  g  is  the  gravitational  acceleration 

Expansion  of  Equation  A. 19  yields. 


P  P^§^(e  +j(u^  +v^))       =      -V-(-pk^VT^)     -V-(pPV) 


-    y    y    ^ p    T.  .  V.    +    hZ  (T     -T    ) 


1    D         1 


c        a' 


+   p  p    V  g 


(A. 20) 


where 


v(pT-v)     =    I  I 


,       .  p     T.    .   V,         =        V(pT        V       +PT        V 

•■    9x,  ij     j  rr    r      '^    rz    z 


1  :       1 


+    pi       V       +    pT       V    ) 
^    zr   r        "^    zz    z 


3  3 

-r— (pi       U  +  PT      v)     +   Tr-(P^       u+pi      v)  -(A. 21) 

3r   '^   rr        ^   rz  3z    '^    zr        ^    zz 


And  so  Equation  A. 19  becomes. 


rDe  ,  D  ,1  ,  2    2^ 
PPa^Dt-^dt^I^^  +^  ) 


]   =   V-  (pk^VT^)  -^-(pT^^u  +p  T^^v) 


a   a    3r  '^  rr 


rz 


~  37^P^zr^  ^P^zz^^  ■  (VpP-V  +PPV-V) 


+  hZ(T  -T  )  +  p  p  V  g  .      (A. 22) 
c    a     ^  a 
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As  in  the  porous  solid  heat  transfer  equation,  the  time  and 
position  dependent  porosity  appears  inside  the  differential 
Expanding  terms,  the  air  energy  equation  becomes, 

P  P^§f  +  -2^  §^(u^  +v^   =   V-  (p  k^VT^)  -  [VpP-v  +p  P  V-V] 

-  IF^P^rr^^P^rz^) 

-  ll^P^rz^^P^zz^)  ^  hZ(T^-T^: 


+  p  p  V  g  .  (A. 23) 

3. 


Consider  the  momentum  equation  for  the  r-direction, 


R  Momentum 

d     ,  V      1,8,       2,    9,       >    9pP 

^(PP^U)   =   ?(-  37^^PPa^  )  -  3¥^PPa^^^  -  9F 


-  (-|-(prT   )  -  ii  +  l-(p  T   ))        (A. 24 

r  9r  ^    rr      r     3z  '^  rz 


Consider  the  momentum  equation  for  the  z-direction. 


Z  Momentum 


It(PP^v)   =   1(-  l-(rpp^uv)  -  |^(p  p^v2)  -  l^(pP) 


-  (pl^^P-rz^  ^k(P-zz)^  ^  PPa^        ^^''' 
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and  the  continuity  equation, 


Continuity 


3  /    \         o  /    \    /     3u  ,     a.  \         3. 
3^(PPa)   =   -  ^^(PPa^  -  (P  Pa  3^  ^  -^"^  "  ""   "9^ 


3v 
P  Pa  97  (^-26) 


Multiplying  the  continuity  equation  through  by  u,  and  substi- 
tuting this  into  Equation  A. 24,  the  r-momentum  equation 
becomes, 

9u  9u  9u         9pP 

PPa9t      =      -PPa^97  -    P  Pa^    97  ■      9F 

-    (-  |-(p  r  T       ) —  +   |-(PT        ))  (A. 27) 

^r    9r   '^        rr  r  9z    ^     rz 

Multiplying  Equation  A. 26  by  v,  and  substituting  this  into 
Equation  A. 25,  the  z-momentum  equation  becomes. 


PP 


9v  9v  9v    9pP    ,13/       N 

a  Tt   =   -PPa^I7  -  P  Pa^  -97  -  gV  ■  ^F  g^^P^^rz^ 

^  k^P^zz^)  ^  PPa^  ^^-'^^ 


Multiplying  Equation  A. 27  by  u  and  noting  that. 


Du  9u  ,       2  9u  ,         9u  ,,    ^^x 

PP^u—  =   pp^u3^+pp^u   —^pp^uvg^  (A. 29) 
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Equation  A. 27  becomes, 

Du      _  8pP  /I    3     /  s  ^^66       9 

A. 30 


PPa^Dt      =      -^    ar  -    ^^F   a^^P  ^^rr^    -   -T-^3^(P  Vz^^ 


Multiplying  Equation  A. 28  by  v  and  noting  that. 


Dv  3v  8v  2    dv  ,       .. . 

PPa^Dt      =     PPa^    3t    ^   P  Pa^   ^    37  ^   PPa^      Jz  ^^'^^^ 


Equation  A. 2  8  becomes 


PPa^  Dt   =   -  ^  3r  -  ^^F  T?(P  ^^rz^)  ^  Sl^P^zz^^  +  PP^^  9 

(A. 32) 


Expanding  Equation  A. 30  yields, 


Du         3pP      ^P^rr    u      ^  ^P'ee      3  , 
P  Pa^  Dt   =   -  ^  3F  -  ^  -JT-  -   FP^rr  ^  -T"  "  ^  3¥^P^rz^ 

(A. 33) 


Expanding  Equation  A. 32  yields, 

Dv                     3pP              ^^P^rz^  V  „    9     ,  ,    ^ 

PPa^Dt      =      -^     aT  -    ^  — 3F FP^rz    -    ^    Sl^P^zz^    ^    PPa^5 

(A. 34) 


The  energy  Equation  A. 23,  after  substituting  and  expanding 
terms,  is, 
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+      U      r +      pi  ^r— 

8r  '^    rz    3r 

^   ^    al^P^zr^    ^   P^zz    97 


^   ^  ll^P^zz^)    -     (^PP-^ 


+    pPV-V)     +PP     vg 

Si 


+    hZ(T     -  T    )  (A. 35) 

c         a 


Substituting  Equation  A. 33  and  Equation  A. 34  into  Equation 
A. 35  yields. 


P  P  -KZ      =      V-  (pk  VT  )  -  p  P  Div  V  +  hZ(T  -T  ) 
a  Dt  a   a  c    a 


,      9u  ,       9v  ,        9u  ,       9v, 
^P  ^rrlF  +  P  ^rz  97  +  P  ^zr  Ti  ^  P  ''zz  97^ 

+  — pT    +  —pi  (A. 36) 

r       r  ^  rr    r  '^  rz 


The  viscous  dissipation  terms  in  Expression  A. 36  are 
neglected  because  the  fluid  is  a  gas  flowing  at  low  velocity 
(see  [Ref.  9]).   Therefore,  the  energy  equation  for  the  air 
in  the  porous  medium  is. 
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P  P  7§-      =   V-  (pk^VT  )  -  p  PDiv  V  +  hZ(T^  -T^: 


;a.37 


with  specific  enthalpy  for  a  gas  defined  by. 


h      =   e  +  - 


:a.38) 


De/Dt  can  be  expressed  as, 


De 
Dt 


Dt 


_1  DP  ^  _^^ 
Pa  °t     p2   Dt 


:a.39: 


where  k    is  the  specific  enthalpy,  and  e  is  the  specific 

energy  (internal  plus  kinetic).   Multiplying  the  continuity 

2 
equation  through  by  P/(pp  )  obtains, 

a 


PP 


^l^(PPa)  ^  ^rlr^PPa^^)  ^  I^ll^PPa^^   =   °-  ^^'''^ 


PP, 


pp. 


Expanding  yields. 


— I^P  TT  "■  Pa  # 
PP. 


PP 


PP. 


■2Ug^CpP^)     - 


a  3 


PP. 


2   r   3r 


(ru) 


T:7(PP.) 2  P  P-^^      (A. 41 


PP. 


2  Bz^^^a' 


PP 


a  3z 


In  vector  notation.  Equation  A. 41  becomes. 
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2  3t    pp   at        2     ^^a   pp 

•"^a  a 


-^  Div  V  (A. 42 

Pa 


Employing  Stokes  (substantial)  derivative  notation,  Equation 
A. 42  becomes, 

P  ^^^  P  Dp     P  ^.   r>,  ,^A■,^ 

-^   -=-—   =   -  — r—  -=^  -    —  Div  V  A. 43) 

7^  Dt         PP^  Dt    p 
Pa             a       a 

Substituting  Equation  A. 43  into  Equation  A. 39  yields. 


D|   =   C|.-AD|-XDivV--^5E  (A. 44 

Dt      Dt    p^  Dt    p  pp   Dt 

a       a  a 


Substituting  Equation  A. 44  into  Equation  A. 37  yields, 


P^a  5^-  P5I-  P  §1   =   VMpk^VT^)  *  hZ(T^-T^)         (A.45) 


Simplification  yields. 


PPail-^   =   ?-(pk^VT^)  +  hZ(T^-T^)  (A.46) 


Invoking  the  Maxwell  relations  for  a  simple  compressible 
substance. 


d/i   =   T  ds  +  -  dP  (A. 47) 

P 
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dT    S 
ds   =   C   ^  -  -  dP  (A. 48) 

p  T    p  ^     ' 


and  recalling  that. 


B   --^|f)p  (A.49) 


the  equation  of  state  for  a  perfect  gas, 


P   =   p  RT  (A. 50 


simplifies  Equation  A.49  such  that, 


RT  ,   -P  >  1  /TV    CI  X 

-Tr( 2^   =   T  (A. 51) 

^  RT        ^ 


Thus,  Equation  A. 4  8  becomes. 


dT    1 

ds   =   C   ^ ^  dP  (A. 52) 

p  T    pT 


Substituting  Equation  A. 52  into  Equation  A. 47  and  cancelling 
terms  yields. 


dh   =   C   dT  (A. 53) 

P 


or 


TM  DT 

R^-   =   C   ^  (A.  54) 

Dt       a  Dt 
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Substituting  Equation  A. 54  into  Equation  A. 45  yields, 

3T 
PPa^a  -^     =   VMpk^VT^)  -  PP^CJV-V)T^ 

+  hZ(T^  -T^)  +  ^^§P-  (A. 55) 

Making  the  assumption  that  pP  changes  very  little  with  time 
[Ref .  9] ,  (this  assumption  was  subsequently  confirmed  by  the 
model) ,  i.e.. 


DE£   =   (v.V)pP   =   u  ^i£^  +  V  li£PI  (A.56) 

Dt  ^  or  9z 


The  final  air  heat  transfer  (energy  conservation)  equation 
is. 


PPa-^aTF  =   '-(Pi^a'Ta*  -  PPa=a'^-'>Ta^  "  lf^+  ^  If^ 

+  hZ(T  -T  )  -  (A. 57) 

c    a 


In  vector  form.  Equation  A. 57  becomes. 


+  (V-V)pP  +  hZ(T  -T  )  (A. 58) 

c    a 


The  expressions  used  to  obtain  the  properties  and  parameters 
in  the  coefficients  of  Equation  A. 57  are  presented  in 
Section  III.E. 
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D.   OXYGEN  MOLECULE  DIFFUSION  EQUATION 

The  final  consideration  in  formulating  the  field  equations 
for  the  model  is  the  transport  of  oxygen  molecules.   The 
oxygen  molecule  transport  equation  is  obtained  by  a  conser- 
vation of  species  balance  on  the  differential  volume  of  air, 
dV  =  pdV.   The  convention  used  for  the  species  balance  into 
a  differential  volume  is, 


Oy    in   =   O2  out  +  0„  consumed  +  0^  accumulated  (A. 59 


The  transport  mechanisms  considered  were  diffusion  due  to 
concentration  gradients  (Fick's  law),  convection,  and  air 
consumption  by  combustion.  The  species  balance  on  oxygen 
becomes. 


pm, . ^^dA    +  pm,.-_dA    +  pm     dA    +  pm     dA 
^   diff         diff         conv       '^    conv 
r  z  r 


pm^.ffdA 


r+dr 


+  pm^.ffdA 


z+dz 


+  pm     dA 
conv 


r+dr 


+  pm     dA 
^  conv 


z+dz 


+  m     dA'  +  pm   dV 
cons         ace 


:a.60) 


Representing  terms  on  the  right  side  by  Taylor  series  expan- 
sions (neglecting  higher  order  terms),  i.e., 


pm,  dA 


£ . +d£ . 
1    1 


3 


=   pm^^dA     +  y^(pmj^dA)d£^ 


(A. 61) 


^i 
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where, 


k  =  conv,   £.  =  r  or  z      dA^  =  r  de dZ  (A. 62) 


and   dA    =   r  dr  d0  (A.  63) 


Then  Equation  A. 61  becomes. 


-^(pmdA^)dr      =      -^(r  p  u  (|)  de  dz)  dr 


=      ^  I— (r  pu  (p)  r  de  dZ    dr  (A. 64) 

r    d  r 


and 


■;r~(p  mdA    )  dz      =      ^^— (p  v  d)  r  dr   d0)dz 

aZ  Z  d  Z 


=      47r(p  v4>)r  de  dz    dr  (A. 65) 

d  Z 


Substituting  Equations  A. 64  and  A. 65  into  Equation  A. 60, 
cancelling  terms  and  rearranging.  Equation  A. 61  becomes. 


-7— -(plmj.^^  +  m     )  dA  )  dr  -  T—{p{m^.^^    +  m     )  dA  )  dz 
9r  ^   diff     conv   r       3z^   diff     conv    z 


-  m     dA'   =   pm    dV  (A. 66) 

cons        '^  ace 


Substituting  the  following  expressions  into  Equation  A, 66, 
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m,  .  ..   =  -    V^  y  (t)  (A. 67 

diff        e   ^ 


m       =  u(t),     m        =   vd>  (A. 68) 

conv,r  conv,z 


m       =   R^  (A. 69 

cons      O^ 


i^      =11-  (A. 70 

ace      8t 


dividing  both  sides  by  dV,  and  letting  dA'/dV  equal  the 
specific  internal  area,  z,    the  oxygen  molecule  diffusion 
equation  becomes, 


V-  (pP^Vcj))  -  V-(P(|)V)  -  Rq  Z   =   £^  (A. 71) 


The  methods  and  expressions  for  obtaining  the  properties  and 
parameters  in  the  coefficients  of  Equation  A. 71  are  presented 
in  Section  lll.G. 
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POROUS    MEDIUM 


AIR 


SOLID 


dV 


0 


+ 


p.dV 


(r-p)-dv 


Figure  A.l 


Separating  A  Differential  Volume  of 
Porous  Medium  into  Respective  Volumes 
of  Solid  and  Air 
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APPENDIX  B 
AUXILIARY  EQUATION  FORMULATION 

A.   CONTINUITY  EQUATION 

The  continuity  equation  for  a  fluid  in  a  porous  medium 
is  expressed  by. 


D(pp^) 

—^ +pp^DivV   =   0  (B.l) 


or  in  equivalent  form, 


3  - 

9^(PP^)  +  (V-V)pp^  +  pp^(V'V)   =  0  (B.2) 


and 


^    <         ^  .    ^^'^  .    ^^^'a^  ^     A  3(ru)  ^  3v, 

^(PP^)  +  u  -^^  +  V  —^^-   +  pp^(-  -^p-  +  ^)  (B.3) 


B.   LAGRANGE  POLYNOMIAL  APPROXIMATIONS  FOR  THERMAL  PROPERTIES 

Relations  for  calculating  the  dynamic  viscosity,  thermal 
conductivity,  and  specific  heat  at  constant  pressure  of  air 
at  different  temperatures  were  required.   Second  order 
Lagrange  polynomial  fits  to  empirical  data  provide  a  simple 
method  for  obtaining  the  relations  required.   Vatikiotis 
[Ref.  9]  gives  the  details  for  arriving  at  the  resulting  set 
of  polynomials. 
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U      =      -3.308  xio"^    T^    +    4.633   x  lo"^    T      +    4.427  x  lo"^  (B.12) 


k        =      -2.608  xio"'^^    T^    +    1.930  x  10~^    T      +    1.361  x  lo"^       (B.13 


C         =      -1.293   X 10    ^    T^    +    2.758   x 10~^    T      +    .238  (B.14) 

a  a  a 


Each  expression  obtains  property  values  within  two  percent 
of  the  data  presented  in  [Ref.  9]  for  temperatures  up  to 
3000  degrees  Fahrenheit. 

C.   CARBON  PARTICLE  SURFACE  RECESSION 

The  following  analysis  of  particle  diameter  consumption 
assumes  that  the  fuel  particle  surface  recedes  uniformly. 
This  assumption  is  reasonable  if  the  particle  diameter  is 
small  in  comparison  to  the  cylinder  geometry,  i.e.,  negligi- 
ble boundary  effects.   In  addition,  since  the  velocities  are 
low,  the  hydrodynamic  effects  on  the  unifoirmity  of  the 
particle  surface  recession  are  negligible.   The  analysis  also 
assumes  there  are  no  significant  thermal  gradients  within 
the  particle,  i.e.,  an  isothermal  carbon  particle.   This  is 
also  reasonable  for  the  small  particles  examined  in  this 
analysis.   A  mass  balance  must  be  performed  on  the  particle 
and  equated  with  the  reaction  rate.   This  equivalence  yields 
the  following  expression, 

^   =   R  Z  (B.15) 

dt       c 

j 
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or  equivalently , 

d(l-p)p 


dt 


Simplification  of  Equation  B.15a  yields, 


-   =   R  2  (B.15a: 


R  2 
|E  =  ^  (B.16) 


Substituting  the  porosity  expression  for  spherical  particles 
yields , 


at<^'5''>    =   ¥  <^-"» 


Application  of  the  time  operator  yields, 


_  3   d^  4   =   _^  (B.li 

^      D^       Pc 


Isolating  the  diameter  time  derivative  yields, 

-  2  R  ZD^ 

d   =   "=—2 (B.19) 

p   IT  d 
c 
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APPENDIX  C 
GALERKIN  FEM  FORMULATION 

A.   FINITE  ELEMENT  METHOD 

The  solution  of  the  system  of  coupled,  nonlinear  partial 
differential  equations  given  by  Equations  3.9,  3.22,  3.26 
and  3.27  —  subject  to  boundary  and  initial  conditions,  was 
obtained  by  a  Galerkin  formulation  of  the  Finite  Element 
Method. 

1 .   Galerkin  Formulation 

A  Galerkin  formulation  of  the  Finite  Element  Method 

was  used  to  obtain  solutions  of  the  porous  solid  and  air 

energy  equations,  the  oxygen  diffusion  equation,  and  the 

continuity  (pressure — Darcy ' s  law)  equation.   A  convenient 

form  of  Equations  3.9,  3.22,  3.26  and  3.27  was  used  in  the 

formulation  where  the  spacial  coordinates,  r  and  z  were 

nondimensionalized  by  £  =  z/z  and  n  =  r/r  . 

o  o 

The  closed  domain  defined  in  (r,z)  space  by  (0,0), 
(0,1),  (1,1),  and  (1,0)  was  partitioned  into  NEL 
(2* (NRNP-1) * (NZNP-1) )  contiguous  area  elements  obtained 
from  a  NSNP  model.   NSNP,  NRNP  and  NZNP  are  the  number  of 
system  nodal  points,  radial  points,  and  axial  points, 
respectively.   The  four  field  variables  T^ ,  T  ,  P  and 
were  approximated  by. 
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NSNP 


T 


=  ipT.(n,£,t)   =    I      N  (n,£)e-.(t)         (c.i 


NSNP 
T^   =  ^'.■ir],e,t)       =  I       N  (n,e)e-.(t)  {C.2: 


j  =  l 


NSNP 


P   =  ilJ^.{r],e,t)       =  I       N  (n,£)e3.(t)  (C.3: 


NSNP 
<t>      =  ii..{r],e,t)       =    I       N  (n,  e)  e.  .  (t)  (C.4 

-'  j  =  l   J       ^J 


where  N.  for  j  =  1,...,NSNP  is  a  set  of  specified  linear  basis 

functions  with  local  support,  and  the  sets  9-,.,  9^.,  9-,., 

^^         '  1] '   2  J    3^ 

and  9..;  j  =  1,...,NSNP,  are  the  solution  coefficients  to  be 

determined.   The  N.  were  selected  to  satisfy  the  condition 

N.(NP.)  =6.-  where  the  Kronecker  delta,  6..,  is  defined  by 
J    1      13  ij  -^ 

6  .  .  =  1  for  i  =  j  ,  and  6  .  .  =  0  for  i  5^  j  .   As  a  result,  9 ,  .  , 

9^.,  9-,.  and  9..  are  the  values  \b^  ,    \b  ^ ,    i) -,    and  \b  .    at  the  nodal 
2 j    3]       4j  r2_/  r 2 '  ^3      ^4 

points  (i.e.,  il) .  ■  (n  /  e  ,  t)  =  9  .  .  (t)  )  . 

Area  interpolation  functions  (shown  in  Figure  C.I) 
were  used  as  the  linear  basis  functions  which  provide  the 
necessary  function  continuity.   As  a  measure  of  error,  a 
residual,  R.,  is  defined  for  each  field  equation  by. 


R.       =       ip    +    A.(i|;)  -  F  +  C[T_(})]*  (C.5) 

1      ~     1       -„    ~,   c 
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where  A.  denotes  the  spatial  operator  of  the  ith  equation 
and  the  asterisk  denotes  that  this  term  appears  in  the  carbon 
temperature  and  oxygen  concentration  equations  only.   The 
term  arises  from  the  reaction  terms  and  is  developed  in 
subsection  3  of  this  section.   The  following  notational 
convention  for  differentiation  is  adopted, 


(   )i   =  ^j^  (C.6) 


(  ■  )   =  ^^  (C.7) 


For  field  equations  3.9,  3.22,  3.26  and  3.21,    the  residuals 
are, 

R^    =   (l-P)PcCc^c  "  V-  (d-p)  (k^)  VT^)  +  hZ(T^-T^) 
c 

-  R  Z  (C.8) 

g 


V  =      PPa^a^a  "  ^*  (P^^a^^a^  "  hZ(T^-Tj 


"A 


+  p  p  C  (V-V)T   -  (V-V)pP  (C.9) 

a  a      a 


R    =   (pp  )  -"V-(— VP)  -  — (VP-V)pp 
p       '^^^a        py   '    py      '  ^^i 


(CIO) 


^^      =  p(^    -    y-  ip  V      V(1))+V-(P({)V)+RqZ  (C.ll) 
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where  the  coefficients  of  the  response  variables  are  them- 
selves functions  of  the  response  variables,  and  thus,  the 
equations  are  nonlinear.   In  accordance  with  the  Galerkin 
method,  the  final  system  of  ordinary  differential  equations 
was  obtained  by  setting  each  residual,  R.,  orthogonal  to 
each  basis  function,  N.,  that  is. 


[NR.  dA=0  (C.12) 

A^  -       "- 

e 


The  4*NSNP  ordinary  differential  equations  given  by  Equations 
C.12  retain  the  character  of  the  original  set  of  partial 
differential  equations,  i.e.,  self-adjoint  operators  yield 
symmetric  matrices  and  non  self-adjoint  operators  yield 
nonsymmetric  matrices.   Thus,  linear  field  operators  trans- 
form to  matrix  operators  and  nonlinear,  coupled  algebraic 
operators.   Incorporation  of  the  boundary  conditions  resulted 
in  4*NSNP  nonlinear  coupled  ordinary  differential  equations. 


F(i,ip,t)   =   Ai|;  +  Bi|;-F  +  C^j[^^  ip  ^]  (C.13) 


subject  to  initial  conditions,  where  A  is  a  (4*NSNP) * (4*NSNP) 
matrix,  B  is  the  matrix  associated  with  the  linear  field 

operators  in  Expression  C.5,  F  is  an  excitation  vector, 

* 
C.  .  is  a  3x9  matrix  arising  from  a  bilinear  operator  treat- 
ment of  the  reaction  terms  in  Expressions  C.8  and  C.ll. 
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Adopting  the  convention, 


NSNP 
^  .       =   N  9   =     y   N.0 .  .  ,   i  =  1,4  (C.14) 


and  performing  an  integration  by  parts  on  the  second  order 
derivatives  yields. 


//nR     dA=-^      (l-p)N(k    ) (T    )     d£    +    //     (1-p) (k    )N     (T    )     dA 
A~c  JJe  ~ecn  ^  e-^rcr 

e  e 


-  //  -^^  N(k JdA  +  //  (1-p)  (k  )N  (T  )  dA 
e  e 


■<j)  (1-p)  (k  )N(T  )  dA  +  //  hZN(T   -T  )  dA 

ie  e  .   c  n     j^^    ^   c    a 

e 


//  R   ZNdA  +  //  p  p   N  T   dA  (C.15) 

A   g  -    A     c  ~  c 
e  e 
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//  NHp  dA     =^^pk    N(T    )     d£    +    //    pk    N     (T    )     dA 
A     ^     a  ci~       d    XI  .  ^  ~i       d    i 

e  e 


A  A 

e  e 


(6         pk    N(T    )     dZ    -    j  (    hZN(T     -T    )  dA 

e 


//    upN(P)     dA    -    f/    u    ^  NPdA 

A  -         ^  A  ^^    ~ 

e  e 


jj    vpN(P)^dA    -    //    V   If  NPdA 

A  ^  A  ~' 

e  e 


+    //    p  p     C     uN(T    )     dA    +    //    p  p     C      vN(T    )     dA 
e  e 


+    //    p  p     C    N  T      dA  (C.16; 

A  a    a~     a 

e 
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p  m 
(In  R  dA=-6       -^(P)  d^  +  //  -^-   N  (P)  dA 


p  m 
-  //  -A_N(p)  dA  -  ^    Pa"^  ,,,,,  ^, 


^  //  !:£ 


P^m  8pp 


N  (P)  dA  +  //  —  -^  N(P)  dA 
A    y   ~z    z      A   py   9^   ~    r 


^  (/  F^^^'P'z^'^  ^  //   R^?P<5A  (C.17 

A   ^  A      a 


//NRdA  =  -  ^   pP  N(4))^d£  +  //  pP  N  (4))  dA 


-    jl    —^   N((}))  dA  -  ^   pP  N((l))  dil  +  //  pP  N  ((}))  dA 

A    ^   ^    ^        ile   ®~    ^      A     ^~^    ^ 
e  e 


(rpu) 
+  //  puN(<|))  dA  +  //  pvN(4))  dA  +  // N({)dA 

A^    ~    "^      K         ^         ^  A      ""    ^ 

e  e  e 


+  //  (pv)  N(})dA  +  //  R   ZNdA  +  //  pN(|)dA  (C.18) 

A       ^~       A    "2  ~      A    ~ 
e  e  e 


The  line  integral  terms  in  each  expression  above  are  boundary 
terms  which  permit  incorporation  of  natural  boundary  conditions 


[ 
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Implementation  of  the  boundary  conditions  is  presented  in 
Section  2  of  this  appendix.   The  coefficients  in  Equations 
C.15,  C.16,  C.17  and  C.18  are  temperature  dependent  properties, 
and  were  taken  as  the  average  values  of  the  properties  over 
an  element.   In  the  limit,  as  the  elements  get  smaller  (i.e., 
NSNP  -»-  «>)  ,  the  average  values  of  the  coefficients  converge 
to  the  exact  values.   Inspection  of  expressions  C.15,  C.16, 
C.17  and  C.18  yields  the  six  operators, 


//  N(i|;)  dA  (C.19) 

A   ~    ^ 
e 


//  Ndi;)  dA  (C.20: 

A   -    2 
e 


//  N  (^)  dA  (C.21 

A   -^   ^ 
e 


//  N  {^)    dA  (C.22 

A   -2    z 
e 


//   NipdA  and  //   Ni|;dA  (C.23 

A      -  A      ~ 

e  e 


N((J;)     dl  (C.24) 

£e~         n 
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To  formulate  these  operators,  the  global  linear 
shape  function,  N.,  are  defined  on  the  local  level  by, 


T        T 

N  e  ^  ?  e  (C.25) 


where  the  natural  coordinates  E,^  ,  E,^    and  E,-,    are  defined  by 
Figure  C.2, 


A. 
^i      =      ^   '.       i  =  1'2,3  (C.26) 

e 


and  A   is  the  area  of  the  elenmet  and  the  A.  are  the  areas 
e  1 

(in  Figure  C.2)  subtended  by  lines  from  a  point  P(r,z) 
inside  the  triangle  to  the  triangle's  vertices  (r-;,z-:, 
j  =  1,2,3).   The  local  shape  functions  (i.e.,  the  elements) 
have  the  following  properties, 

^ij   "   ^ij'^-    '   ^  "  ■^'^'      ^    "  ^'^  (C.27) 

5  .(NP  )   =   6. .  ;    i  =  1,3,   j  =  1,3  (C.28) 

J        -'-  J 

Having  defined  the  local  shape  functions,  the  ele- 
mental matrix  operations  C.19  through  C.24  are, 

//    N(ip)        dA     =       |tl^i]?^lt  ^^-25^ 

A       ~ 
e 

//    Ndp)       dA     =      It^.LQelt  ^^'^^^ 


A 
e 
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//  N^(i|i)  dA 

A       ^ 
e 


4A 


b.b.] 9  ,^ 
■  1  J  ~elt 


(C.31) 


//  N  (^)^dA  =   [c  c  J 


A 


z  '  z 


4A 


i  j  ~elt 


:c.32 


A   " 
®  k  =  2,  i  =  j 


(C.33) 


<j)         kN(i|;)  d£   =    ^    kN(tJ;)  dr  -  ^    kN(4;)  dz 


:c.34) 


where  a(i,j)  coefficients  of  the  element  matrices  A(3  x3)  are 
given  by, 


a(i,j)   =   [g]  ;     g=g(i,j) 


:c.35) 


and  the  vector  9  -,   is  given  by, 


9   TV 

~elt 


(C.36: 


The  derivations  of  these  operators  are  presented  in  Section 
A.  5  of  this  appendix. 


116 


2 .   Implemen-tation  of  Boundary  Conditions 

Having  formulated  the  system  matrices  for  the  field 
equations,  treatment  of  the  boundary  conditions  is  now  dis-  . 
cussed.   Each  field  equation  is  treated  individually. 

a.   Porous  Solid  Energy  (Heat  Transfer)  Equation 

There  are  three  types  of  boundary  conditions  that 
can  occur:   Dirichlet,  Neumann  and  Cauchy  (mixed)  boundary 
conditions.   The  treatment  of  each  is  as  follows.   For 
Dirichlet  boundary  conditions,  one  may  specify  an  equation 
to  be  a  linear  equation  (i.e.,  independent  of  time)  by  placing 
a  1  on  the  diagonal  of  the  system  stiffness  matrix  correspond- 
ing to  the  particular  degree  of  freedom  at  hand  and  setting 
each  time  derivative  matrix  coefficient  for  the  same  equation 
equal  to  zero.   An  alternative  scheme  involves  setting  the 
time  derivative  diagonal  term  equal  to  one  and  setting  all 
stiffness  coefficients  (in  the  same  equation)  equal  to  zero. 
In  this  manner  one  specifies  the  Dirichlet  value  as  an  initial 
condition  whose  residual  is  identically  zero  and  thus  is 
invariant  with  time.   Treatment  of  Neumann  boundary  conditions 
is  as  follows. 


(f>        kVxdC  =      <f>        p  d^   ^   -  -^  added  to  local  F(l) 


ie  le 


2 

and  local  F(2)  at  the  element 
nodal  points  (C.37) 


In  the  above  expression,  x  is  a  response  variable  and  p  is 
an  average  value  of  flux  across  an  element.   Treatment  of      ^ 
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Cauchy  boundary  conditions  is  as  follows, 


le.      1 

2 

1 

^1 

^         kVxd?      = 

=      j>        k(x 

-x^) 

-> 

ie 

£e 

•~  V--               J 

1 

2 

X 
^  2 

—    le. 
scattered  into  3   and   -  K,   ^r  added  to 

/e  2 


F(l)  and  F(2)  at  the  element  nodal  points 


;c.3i 


For  a  complete  development  of  the  theory  for  incorporating 
boundary  conditions,  see  [Ref.  28]. 

3 .   Treatment  of  the  Reaction  Rate  Term 

An  exponential  reaction  rate  term  appears  in  both 
the  porous  solid  and  oxygen  diffusion  equation.   The  reaction 
expression  is, 


R^   =   A  4)''exp(-E/R^T^) 


(C.39) 


In  the  carbon  equation.  Expression  C.39  appears  as. 


R  Z   =   (R  AH  )Z 


(C.40) 


In  the  oxygen  concentration  equation.  Expression  C.39 
appears  as. 


^2 


(C.41) 
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The  reaction  rate  term  may  be  treated  in  various  ways. 

a.  Treatment  as  an  Excitation  (Force  Term) 

This  treatment  is  the  easiest  of  all  schemes. 
The  term  is  merely  evaluated  at  the  last  time  step  and  is 
assumed  to  be  constant  over  the  next  integration  time  step, 
i.e., 

R^   =   {A  R^  exp(E/RT^)}  (C.42) 

The  superscript  *  indicates  evaluation  occurring  at  the 
previous  time  step.   The  value  of  the  term  is  evaluated 
at  the  ith  nodal  point  and  inserted  as  a  factor  in  the 
corresponding  carbon  energy  or  the  oxygen  diffusion  equation 
of  the  ith  nodal  point. 

b.  Linear  Operator  Treatment  of  O2  Concentration 

In  order  to  realize  an  improvement  over  the  first 
treatment,  one  may  retain  a  portion  of  the  reaction  expression 
as  an  operator  by  making  the  following  rearrangement: 


R^   =   (A  R^  (t)""^exp(-E/R  T^)  }*•  [<})]  (C.43) 

c  0*5  u  c 


where  [(p]    denotes  a  spatial  operator  treatment  of  the  response 

variable. 

c.   Temperature  and  0^  Concentration  Bilinear 

Operator  Treatment  of  the  Reaction  Rate  Term 

The  bilinear  operator  treatment  of  the  reaction 

rate  is  the  present  method  of  treatment.   The  reaction  rate 
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term  is  rearranged  as  follows, 


R^   =    ^ ^ [T^$]       (C.44 


Letting 


T^   =   §^e^  (C.45) 


and 


(|)   =   §  e^  (C.46) 


and  invoking  natural  coordinates  [Ref.  28],  an  elementally 
averaged  contribution  results  in  the  following  area  integral, 


c  =  c  //  5  5  e,  ?  e,  dA 

e 


(C.47: 


Expanding  yields, 

C   =   C  //  l*^   <   ?i5j  {  ?2^j  |^3?j   >   ^^li'^4j^  '       ^^-^^^ 

i  =  1,3,   j  =  1,3 

A  final  explicit  expansion  of  the  integral  looks  like. 
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c  =  c 


3     2      2       2        2 


2       2 

^1^2  ^1^2 


^1^2^3 


^1^2 


2  2 

^1^3  ^1^2^3  ^1^3 


;  ^1^2^3 


.3 
^2 


hh 


2^3  ^2^3 


^Ih 


^1^3 


?l52?3  «1?3 


2        2 

^1^2^3  ^2^3    ^2^3 


^2^3 


^ 


3 
(C.49) 


J 


Invoking  the  integral  formula  of  natural  coordinates  [Ref. 
28]  , 


//    ?i    ^2    ?3    ^^ 
A 


a!b!c!2A 


(a  +b  +c  +  2) 


;c.50) 


results    in. 


C      = 


CA 
e 

60 


6      2      2     J    2      2      1 
2      2      1*262 


2  12 
12  2 
2      2      6 


^Ict 


(C.51) 


where. 


-c     -[ 


(AH 


R 


or 


^H^' 


(C.52) 


The  asterisk  denotes  as  applicable  for  T   or  cj)  equations. 

4 .   Implementation  of  Reaction  Term  in  the  Numerical 
Method 

Franke ' s  (modified  Gear)  integration  routine  requires 

a  calculation  of  (or  an  approximation  to)  the  Jacobian  matrix 
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from  the  system  matrices,  A(t),  B(t)  and  C(t).   In  order 
to  include  effects  of  T   and  (p    arising  from  reaction  rate 
terms  in  the  Jacobian  vector,  and  thereby  improve  the  effi- 
ciency of  the  integration  routine,  the  combustion  terms  are 
incorporated  into  the  residual  equations  (in  DIFMOD) . 
Modifications  to  the  Jacobian  matrix  are  accomplished  in  the 
JACMOD  and  NUITSL  routines.   Reaction  rate  terms  of  Expres- 
sions C.40  and  C.41  contribute  nine  terms  to  the  respective 
residual  equations  and  twenty-seven  terms  to  the  Jacobian 
vector.   Reaction  terms  are  generally  expressed  as, 


3    3 

RT      =         I         I      C^^iii^^^^^.)     ,         k  =  1,9  (C.53 

i=l  j  =  l  -^ 


The  Jacobian  is  defined  as, 

J       E         m±i±LlL   ^   ^C^,^,t)  (C.54) 

Contributions    arising    from   RT    as    a    result   of    9F/3i];   are. 


It^      =         I      C.,  *    ^..  (C.55: 

^li  ;]  =  1  -^ 


where  the  k*  are  compatible  with  the  column  locations  in  the 
C  matrix  of  the  ip,  .  products,  and 
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3  RT 
^4j      1=1 


where  the  k*  are  compatible  with  the  column  locations  in  the 
C  matrix  of  the  ip..    products.   Terms  from  Equation  C.55 
may  be  incorporated  into  the  PW  Jacobian  vector  (containing 
contributions  to  the  Jacobian  from  the  linear  spatial 
operators)  as  contributions  from  Equation  C.8  with  combustion. 
Similarly,  terms  from  Equation  C.56  may  be  incorporated  into 
PW  for  terms  arising  from  the  0^  residual  equations  with 
combustion.   The  remaining  contributions  to  the  Jacobian 
vector,  Equation  C.56  for  the  carbon  equations  and  terms 
arising  from  Equation  C.55  for  the  0^  equations  are  stored 
in  the  PWMOD  (NDOF/2,9)  matrix.   The  column  number  or  variable 
of  differentiation  array,  INMOD  (NDOF/2,9)  and  PWMOD  array 
are  communicated  to  the  NUITSL  routine  so  additional  Jacobian 
terms  not  assimilated  into  the  PW  array  (by  virtue  of  the 
storage  scheme  selected)  may  be  taken  into  account  during 
the  convergence  sequence  for  the  new  iterates.   The  iterating 
scheme  is  of  Newton-Raphson  type.   Thus,  the  final  synthesis 
of  the  Jacobian  including  effects  of  all  terms  arising  from 
combustion  is  consummated. 

5.   Derivation  of  the  FEM  Operators 

In  the  section  on  the  finite  element  formulation, 
the  following  six  differential  operators  were  identified. 
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//  N(-i;)  dA  (C.57) 

A^  ~    ^ 


//  N(ij;)  dA  (C.58) 

A       -  ^ 


e 


//  N  (i^)  dA  (C.59) 


//  N  (.1^)  dA  (C.60) 

A   -2    z 


e 


//  N(ij;)  dA     and     //  N  ip  dA  (C.61) 


A   ^  A 

e  e 


j)        kN(i]j)  di  (C.62) 


where  N.  are  the  global  basis  functions.   These  operators 
are  constructed  on  the  element  level  by  introducing  the 
corresponding  element  basis  functions,  E, .  .      The  global  and 
element  basis  functions  are  related  by. 


\l}.       =      n'^O   ^   ?'^e  (C.63) 


The  derivation  of  the  local  elemental  matrices  (using 
the  local  coordinate  system  depicted  in  Figure  C.3)  according 
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to  the  Galerkin  method  for  the  global  operations  proceeds 
as  follows: 

For  Operator  C.29(also  C.57), 


Global 


Local 


//  ^{ip)  dA 

A   ^ 


//  ?  re  dA 

A   ~  -r- 


:c.64) 


Noting  that, 


3? 


9r 


3       - 


b.0  . 
2A 


and 


9^^ 


9z 


1   = 


c  .9  . 

2A   ' 
e 


J  =  1.3 


(C.65) 


where  the  repeated  index  implies  Einsteinian  notation,  the 
elemental  matrix  becomes, 


A      e" 


(C.66) 


Expanding, 


// 


■^il 


Ll3 


b  .  T 
e 


P 


^   ^2   ^3 


^   ^2   ^3 
^1   ^2   ^3 


e    (C.67) 


For  Operator  C.30  (also  C.58), 


Global 


Local 


//  N{^)  ^  dA 


//  C  5'edA 

A   ~  "^^ 

e 


(C.68) 
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Substituting  the  local  shape  functions  gives, 


C  .  T 

//  ?  (2^)  e  dA 

A   ^    e 
e    ~~~ 


;c.69 


the  elemental  matrix  becomes. 


// 


^1 

^2 

s. 

C  .  T 

^)    e  dA 


2A 


c  c  c 

1    2    3 

c  c  c 

^1  ^2    3 

c  c  c 

^1    2    3^ 


:c.7o 


For  Operator  C.31(also  C.59) 


Global 


Local 


//  N  i^p)  ^    dA 


A 


r    r    e 


//  5^?^e  dA 

e 


:c.7i) 


Substituting  the  local  shape  functions  gives. 


// 


b.  b. 

_X  1 

2A  2A 

e  e 


^^1     ^^2   ^1^3 


dA 


// 


4Af  A 
e   e 


b^b^   b2 
^V3   ^3^ 


^2^3 


edA    (C.72) 


the  elemental  matrix  becomes. 


//  N   (4^)   dA 
A   -r     r 

e 


4A 


b^b^ 
^1^3 


^1^2 


^2^3 


^1^3 

^2^3 

b.   , 


9      (C.73) 
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For  Operator  C.32  (also  C.60), 


Global 


Local 


//   N  di;)  dA 

A    -^    ^ 
e 


//  ?2  ?2  ®  ^^ 


(C.74) 


Substituting  the  local  shape  functions  gives, 


// 

A 


^   ^r-P-  6  dA   = 


2A    2A 
e     e 

e  ~-~~      ~~-~ 


// 


4A    A 
e    e 


C   C        C  C 
^12   ^13 


^iS   S 


C2C3 


^1^3   ^2^3   ^3 


0 dA   (C.75) 


the  elemental  matrix  becomes. 


//  ?2  ?^  ®  "^^ 
e 


4A 


^iS    ^iS 


^1^2   ^2 


C2C3 


"^1*^3   ^2^3   ^3 


(C.76) 


For  Operator  C. 33  (also  C.61) 


Global 


Local 


/  /  N  ij;  dA 
A   ~ 


//  5  ?   e  dA 


;C.77) 


Substituting  the  local  shape  functions  gives. 
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A 

e 


e  dA    = 


// 

A 


^1  hh 

^1^3       ^2^3 


^1^3 
2 


? 


J 


e    dA 


:C.7: 


the    elemental   matrix   becomes, 


//    C    ?      e    dA 
A      ^    ^       ~ 


A 
e 

12 


2  11 
12  1 
112 


(C.79) 


The  last  operator  C.34  (also  C.62)  is  incorporated  into  the 
excitation  vector  as  described  in  the  FEM  formulation  and 
has  been  addressed  in  Subsection  2,  Implementation  of 
Boundary  Conditions. 
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Figure  C.l    Area  Interpolation  Functions 
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-*-  r 


Figure  C.2    Area  Coordinates 
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k 


Figure  C.3    Local  Coordinates 
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